Tim Xiao
email: tim.xiao@cibc.com

An Accurate Solution for Credit Value Adjustment (CVA) and Wrong Way Risk

This paper presents a Least Square Monte Carlo approach for accurately calculating credit value adjustment (CVA). In contrast to previous studies, the model relies on the probability distribution of a default time/jump rather than the default time itself, as the default time is usually inaccessible.

As such, the model can achieve a high order of accuracy with a relatively easy implementation.

We find that the valuation of a defaultable derivative is normally determined via backward induction when their payoffs could be positive or negative. Moreover, the model can naturally capture wrong or right way risk.

For years, a widespread practice in the industry has been to mark derivative portfolios to market without taking counterparty risk into account. All cash flows are discounted using the LIBOR curve. But the real parties, in many cases, happen to be of lower credit quality than the hypothetical LIBOR party and have a chance of default.

As a consequence, the International Accounting Standard (IAS) 39 requires banks to provide a fair-value adjustment due to counterparty risk. Although credit value adjustment (CVA) became mandatory in 2000, it received a little attention until the recent financial crises in which the profit and loss (P&L) swings due to CVA changes were measured in billons of dollars.

Interest in CVA began to grow. Now CVA has become the first line of defense and the central part of counterparty risk management.

CVA not only allows institutions to move beyond the traditional control mindset of credit risk limits and to quantify counterparty risk as a single measurable P&L number, but also offers an opportunity for banks to dynamically manage, price and hedge counterparty risk. The benefits of CVA are widely acknowledged. Many banks have set up internal credit risk trading desks to manage counterparty risk on derivatives.

The earlier works on CVA are mainly focused on unilateral CVA that assumes that only one counterparty is defaultable and the other one is default-free. The unilateral treatment neglects the fact that both counterparties may default, i.e., counterparty risk can be bilateral. A trend that has become increasingly relevant and popular has been to consider the bilateral nature of counterparty credit risk. Although most institutions view bilateral considerations as important in order to agree on new transactions, [START_REF] Hull | CVA and wrong way risk[END_REF] argue that bilateral CVA is more controversial than unilateral CVA as the possibility that a dealer might default is in theory a benefit to the dealer. CVA, by definition, is the difference between the risk-free portfolio value and the true (or risky or defaultable) portfolio value that takes into account the possibility of a counterparty's default. The risk-free portfolio value is what brokers quote or what trading systems or models normally report. The risky portfolio value, however, is a relatively less explored and less transparent area, which is the main challenge and core theme for CVA. In other words, central to CVA is risky valuation.

In general, risky valuation can be classified into two categories: the default time approach (DTA) and the default probability approach (DPA). The DTA involves the default time explicitly. Most CVA models in the literature [START_REF] Brigo | Bilateral counterparty risk valuation with stochastic dynamical models and application to Credit Default Swaps[END_REF], [START_REF] Lipton | Credit value adjustment for credit default swaps via the structural default model[END_REF], Pykhtin and Zhu (2006) and [START_REF] Gregory | Being two-faced over counterparty credit risk[END_REF], etc.) are based on this approach.

Although the DTA is very intuitive, it has the disadvantage that it explicitly involves the default time. We are very unlikely to have complete information about a firm's default point, which is often inaccessible (see [START_REF] Duffie | Swap rates and credit quality[END_REF], Jarrow and Protter (2004), etc.).

Usually, valuation under the DTA is performed via Monte Carlo simulation. On the other hand, however, the DPA relies on the probability distribution of the default time rather than the default time itself. Sometimes the DPA yields simple closed form solutions.

The current popular CVA methodology (Pykhtin and Zhu (2006) and [START_REF] Gregory | Being two-faced over counterparty credit risk[END_REF], etc.) is first derived using DTA and then discretized over a time grid in order to yield a feasible solution. The discretization, however, is inaccurate. In fact, this model has never been rigorously proved. Since CVA is used for financial accounting and pricing, its accuracy is essential.

Moreover, this current model is based on a well-known assumption, in which credit exposure and counterparty's credit quality are independent. Obviously, it can not capture wrong/right way risk properly.

In this paper, we present a framework for risky valuation and CVA. In contrast to previous studies, the model relies on the DPA rather than the DTA. Our study shows that the pricing process of a defaultable contract normally has a backward recursive nature if its payoff could be positive or negative.

An intuitive way of understanding these backward recursive behaviours is that we can think of that any contingent claim embeds two default options. In other words, when entering an OTC derivatives transaction, one party grants the other party an option to default and, at the same time, also receives an option to default itself. In theory, default may occur at any time. Therefore, the default options are American style options that normally require a backward induction valuation.

Wrong way risk occurs when exposure to a counterparty is adversely correlated with the credit quality of that counterparty, while right way risk occurs when exposure to a counterparty is positively correlated with the credit quality of that counterparty. For example, in wrong way risk exposure tends to increase when counterparty credit quality worsens, while in right way risk exposure tends to decrease when counterparty credit quality declines. Wrong/right way risk, as an additional source of risk, is rightly of concern to banks and regulators. Since this new model allows us to incorporate correlated and potentially simultaneous defaults into risky valuation, it can naturally capture wrong/right way risk.

The rest of this paper is organized as follows: Section 2 discusses unilateral risky valuation and unilateral CVA. Section 2 elaborates bilateral risky valuation and bilateral CVA.

Section 3 presents numerical results. The conclusions are given in Section 4. . All proofs and a practical framework that embraces netting agreements, margining agreements and wrong/right way risk are contained in the appendices.

Unilateral Risky Valuation and Unilateral CVA

We consider a filtered probability space (  , F ,   0

 t t F
, P ) satisfying the usual conditions, where  denotes a sample space; F denotes a  -algebra; P denotes a probability measure;   0

 t t F
denotes a filtration.

The default model is based on the reduced-form approach proposed by [START_REF] Duffie | Modeling term structure of defaultable bonds[END_REF] and Jarrow and Turnbell (1994), which does not explain the event of default It is well-known that the survival probability from time t to s in this framework is defined by

            s t du u h Z t s P s t p ) ( exp ) , | ( : ) , (   (2a)
The default probability for the period (t, s) in this framework is defined by

               s t du u h s t p Z t s P s t q ) ( exp 1 ) , ( 1 ) , | ( : ) , (   (2b)
  t F T F X T t D E t V ) , ( ) (  (3a)
where

         du u r T t D T t ) ( exp ) , ( (3b)
where   ).

Next, we turn to risky valuation. In a unilateral credit risk case, we assume that party A is default-free and party B is defaultable. Risky valuation can be generally classified into two categories: the default time approach (DTA) and the default probability (intensity) approach (DPA).

The DTA involves the default time explicitly. If there has been no default before time T (i.e.,

T  
), the value of the contract at T is the payoff T X . If a default happens before T (i.e.,

T t  
), a recovery payoff is made at the default time  as a fraction of the market value2 given by ) (

V

where  is the default recovery rate and

) ( V
is the market value at default.

Under a risk-neutral measure, the value of this defaultable contract is the discounted expectation of all the payoffs and is given by

    t T T T V t D X T t D E t V F | 1 ) ( ) , ( 1 ) , ( ) (          (4)
where Y  is an indicator function that is equal to one if Y is true and zero otherwise.

Although the DTA is very intuitive, it has the disadvantage that it explicitly involves the default time/jump. We are very unlikely to have complete information about a firm's default point, which is often inaccessible. Usually, valuation under the DTA is performed via Monte Carlo simulation.

The DPA relies on the probability distribution of the default time rather than the default time itself. We divide the time period (t, T) into n very small time intervals ( t  ) and assume that a default may occur only at the end of each very small period. In our derivation, we use the approximation

  y y  1 exp
for very small y. The survival and the default probabilities for the period ( t , t t   ) are given by

  t t h t t h t t t p t p          ) ( 1 ) ( exp ) , ( : ) ( ˆ (5a)   t t h t t h t t t q t q          ) ( ) ( exp 1 ) , ( : ) ( ˆ (5b)
The binomial default rule considers only two possible states: default or survival.



. Under a risk-neutral measure, the value of the asset at t is the expectation of all the payoffs discounted at the risk-free rate and is given by

         t t t t V t t y E t t V t q t t p t t r E t V F F ) ( ) ( exp ) ( ) ( ) ( ) ( ) ( exp ) (             (6) where   ) ( ) ( ) ( 1 ) ( ) ( ) ( t c t r t t h t r t y      
denotes the risky rate and

  ) ( 1 ) ( ) ( t t h t c    is called the (short) credit spread.
Similarly, we have

    t t t t V t t t y E t t V            F ) 2 ( ) ( exp ) ( (7)
Note that

  t t y   ) ( exp is t t F   -measurable.
                t i t t t t t t V t t i t y E t t V t t t y E t t y E t t V t t y E t V F F F F ) 2 ( ) ) ( exp ) 2 ( ) ( exp ) ( exp ) ( ) ( exp ) ( 1 0                          (8)
By recursively deriving from t forward over T and taking the limit as t  approaches zero, the risky value of the asset can be expressed as

               t T t T V du u y E t V F ) ( ) ( exp ) ( (9)
We may think of ) (u y as the risk-adjusted short rate. Equation ( 9) is the same as Equation ( 10) in [START_REF] Duffie | Modeling term structure of defaultable bonds[END_REF], which is the market model for pricing risky bonds.

Using the DPA, we obtain a closed-form solution for pricing an asset subject to credit risk. Other good examples of the DPA are the CDS model proposed by J.P. [START_REF] Morgan | The J. P. Morgan guide to credit derivatives, Risk Publications[END_REF] and a more generic risky model presented by Xiao (2013a).

In theory, a default may happen at any time, i.e., a risky contract is continuously defaultable. This Continuous Time Risky Valuation Model is accurate but sometimes complex and expensive. For simplicity, people sometimes prefer the Discrete Time Risky Valuation Model that assumes that a default may only happen at some discrete times. A natural selection is to assume that a default may occur only on the payment dates. Fortunately, the level of accuracy for this discrete approximation is well inside the typical bid-ask spread for most applications (see O'Kane and Turnbull ( 2003)). From now on, we will focus on the discrete setting only, but many of the points we make are equally applicable to the continuous setting.

For a derivative contract, usually its payoff may be either an asset or a liability to each party. Thus, we further relax the assumption and suppose that may be positive or negative.

In the case of , the survival value is equal to the payoff T X and the default payoff is a fraction of the payoff

T X  . Whereas in the case of 0  T X
, the contract value is the payoff itself, because the default risk of party B is irrelevant for unilateral risky valuation in this case.

Therefore, we have

Proposition 1: The unilateral risky value of the single-payment contract in a discrete-time setting is given by

  t F T X T t F E t V ) , ( ) (  (10a)
where

    ) ( 1 ) , ( 1 1 ) , ( ) , ( 0 T T t q T t D T t F T X      (10b)
Proof: See the appendix.

Here

) , ( T t F can be regarded as a risk-adjusted discount factor. Proposition 1 says that the unilateral risky valuation of the single payoff contract has a dependence on the sign of the payoff. If the payoff is positive, the risky value is equal to the risk-free value minus the discounted potential loss. Otherwise, the risky value is equal to the risk-free value.

T X 0  T X
Proposition 1 can be easily extended from one-period to multiple-periods. Suppose that a defaultable contract has m cash flows. Let the m cash flows be represented as 1 X ,…, m X with payment dates 1 T ,…, m T . Each cash flow may be positive or negative. We have the following proposition.

Proposition 2: The unilateral risky value of the multiple-payment contract is given by

           m i t i i j j j X T T F E t V 1 1 0 1 ) , ( ) ( F (11a)
where properly where counterparty credit quality and market prices may be correlated. This type of problem can be best solved by working backwards in time, with the later risky value feeding into the earlier ones, so that the process builds on itself in a recursive fashion, which is referred to as backward induction. The most popular backward induction valuation algorithms are lattice/tree and least square Monte Carlo.

0 T t  and     ) ( 1 ) , ( 1 1 ) , ( ) , ( 1 1 0 )) ( ( 1 1 1 1            j j j T V X j j j j T T T q T T D T T F j j  ( 11b
For an intuitive explanation, we can posit that a defaultable contract under the unilateral credit risk assumption has an embedded default option (see [START_REF] Sorensen | Pricing swap default risk[END_REF]).

In other words, one party entering a defaultable financial transaction actually grants the other party an option to default. If we assume that a default may occur at any time, the default option is an American style option. American options normally have backward recursive natures and require backward induction valuations.

The similarity between American style financial options and American style default options is that both require a backward recursive valuation procedure. The difference between them is in the optimal strategy. The American financial option seeks an optimal value by comparing the exercise value with the continuation value, whereas the American default option seeks an optimal discount factor based on the option value in time.

The unilateral CVA, by definition, can be expressed as

              m i t i i j j j i F X T T F T t D E t V t V t CVA 1 1 0 1 ) , ( ) , ( ) ( ) ( ) ( F (12)
Proposition 2 provides a general form for pricing a unilateral defaultable contract.

Applying it to a particular situation in which we assume that all the payoffs are nonnegative, we derive the following corollary:

Corollary 1: If all the payoffs are nonnegative, the risky value of the multiple-payments contract is given by

           m i t i i j j j X T T F E t V 1 1 0 1 ) , ( ) ( F ( 13a
)
where

0 T t  and     ) ( 1 ) , ( 1 ) , ( ) , ( 1 1 1 1        j j j j j j j T T T q T T D T T F  (13b)
The proof of this corollary is easily obtained according to Proposition 2 by setting

  0 ) ( 1 1     j j T V X
, since the value of the contract at any time is also nonnegative.

The CVA in this case is given by

                   m i t i i j j j j i F X T T T q T t D E t V t V t CVA 1 1 0 1 1 )) ( 1 )( , ( 1 1 ) , ( ) ( ) ( ) ( F  (14)
The current popular CVA model (e.g., equation ( 17) in [START_REF] Pykhtin | A guide to modeling counterparty credit risk[END_REF] and equation (3) in Gregory ( 2009)) is quite different from above either equation ( 12) or equation ( 14).

As a matter of fact, the current CVA model has never been rigorously proved. In order to reflect the economic value of counterparty credit risk, to measure the profit and loss of a bank and to provide proper incentives to traders, a good CVA model must be not only rigorous and accurate but also feasible to implement.

Bilateral Risky Valuation and Bilateral CVA

There is ample evidence that corporate defaults are correlated. The default of a firm's counterparty might affect its own default probability. Thus, default correlation and dependence arise due to the counterparty relations. Default correlation can be positive or negative. The effect of positive correlation is usually called contagion, whereas the latter is referred to as competition effect.

Two counterparties are denoted as A and B. The binomial default rule considers only two possible states: default or survival. Therefore, the default indicator j Y for party j (j=A, B) follows a Bernoulli distribution, which takes value 1 with default probability j q and value 0 with survival probability j p , i.e., The joint probability representations are given by

AB B A B A p p Y Y P p       ) 0 , 0 ( : 00 (15a) AB B A B A q p Y Y P p       ) 1 , 0 ( : 01 (15b) AB B A B A p q Y Y P p       ) 0 , 1 ( : 10 (15c) AB B A B A q q Y Y P p       ) 1 , 1 ( : 11 (15d)
where

j j q Y E  ) ( , j j j q p  2  ,   B B A A AB B A AB B B A A AB p q p q q Y q Y E           ) )( ( :
where AB  denotes the default correlation coefficient and AB  denotes the default covariance. .

State 0 , 0   B A Y Y 0 , 1   B A Y Y 1 , 0   B A Y Y 1 , 1   B A Y Y Comments A & B survive A defaults, B survives A survives, B defaults A & B default
 T X T X T B X  T B X  T AB X  0  T X T X T A X  T A X  T AB X  3
There are two default settlement rules in the market. The one-way payment rule was specified by the early ISDA master agreement. The non-defaulting party is not obligated to compensate the defaulting party if the remaining market value of the instrument is positive for the defaulting party. The two-way payment rule is based on current ISDA documentation. The non-defaulting party will pay the full market value of the instrument to the defaulting party if the contract has positive value to the defaulting party.

Suppose that a financial contract that promises to pay a T X from party B to party A at maturity date T, and nothing before date T where t T  . The payoff T X may be positive or negative, i.e. the contract may be either an asset or a liability to each party. All calculations are from the perspective of party A.

At time T, there are a total of four ( 4 2 2  ) possible states shown in Table 1. The risky value of the contract is the discounted expectation of the payoffs and is given by the following proposition.

Proposition 3: The bilateral risky value of the single-payment contract is given by

      t t F F T A X B X T X T t k T t k T t D E X T t K E t V T T ) , ( 1 ) , ( 1 ) , ( ) , ( ) ( 0 0      (16a) where   ) ( ) ( ) ( 1 ) , ( ) , ( ) , ( ) ( ) , ( ) , ( ) ( ) , ( ) , ( ) ( ) , ( ) , ( ) , (
T T T T t T t q T t q T T t q T t p T T t p T t q T T t p T t p T t k AB B B AB A B AB A B B A B B A B B                (16b)   ) ( ) ( ) ( 1 ) , ( ) , ( ) , ( ) ( ) , ( ) , ( ) ( ) , ( ) , ( ) ( ) , ( ) , ( ) , (
T T T T t T t q T t q T T t q T t p T T t p T t q T T t p T t p T t k AB A A AB A B AB B A A B A A A B A                (16c)
Proof: See the appendix.

We may think of ) , ( T t K as the risk-adjusted discount factor. Proposition 3 tells us that the bilateral risky price of a single-payment contract can be expressed as the present value of the payoff discounted by a risk-adjusted discount factor that has a switching-type dependence on the sign of the payoff.

Using a similar derivation as in Proposition 2, we can easily extend Proposition 3 from one-period to multiple-periods. Suppose that a defaultable contract has m cash flows. Let the m cash flows be represented as . Proposition 4 says that the pricing process of a multiple-payment contract has a backward nature since there is no way of knowing which risk-adjusted discounting rate should be used without knowledge of the future value. Only on the maturity date, the value of the contract and the decision strategy are clear. Therefore, the evaluation must be done in a backward fashion, working from the final payment date towards the present. This type of valuation process is referred to as backward induction.

           m i t i i j j j X T T K E t V 1 1 0 1 ) , ( ) ( F (17a) where 0 T t  and   ) , ( 1 ) , ( 1 ) , ( ) , ( 1 0 )) ( ( 1 0 )) ( ( 1 1 1 1 1 1               j j A T V X j j B T V X j j j j T T k T T k T T D T T K j j j j (17b)
There is a common misconception in the market. Many people believe that the cash flows of a defaultable financial contract can be priced independently and then be summed up to give the final risky price of the contract. We emphasize here that this conclusion is only true of the financial contracts whose payoffs are always positive. In the cases where the promised payoffs could be positive or negative, the valuation requires not only a backward recursive induction procedure, but also a strategic selection of different discount factors according to the market value in time. This coupled valuation process allows us to capture correlation between counterparties and market factors.

The bilateral CVA of the multiple-payment contract can be expressed as

                  m i t i i j j j t i i F X T T K E X T t D E t V t V t CVA 1 1 0 1 ) , ( ) , ( ) ( ) ( ) ( F F (18)

Numerical Results

In this section, we present some numerical results for CVA calculation based on the theory described above. First, we study the impact of margin agreements on CVA. The testing portfolio consists of a number of interest rate and equity derivatives. The number of simulation scenarios (or paths) is 20,000. The time buckets are set weekly. If the computational requirements exceed the system limit, one can reduce both the number of scenarios and the number of time buckets. The time buckets can be designed fine-granularity at the short end (e.g., daily and then weekly) and coarse-granularity at the far end (e.g. monthly and then yearly). The rationale is that the calculation becomes less accurate due to the accumulated error from simulation discretization, and inherited errors from calibration of the underlying models, such as those due to the change of macro-economic climate. The collateral margin period of risk is assumed to be 14 days (2 weeks).

For risk-neutral simulation, we use a Hull , the CVA value actually decreases while the threshold A H increases. This reflects the bilateral impact of the collaterals on the CVA. The impact is mixed in Table 4 when both parties have finite collateral thresholds. Next, we examine the impact of wrong way risk. Wrong way risk occurs when exposure to a counterparty is adversely correlated with the credit quality of that counterparty, while right way risk occurs when exposure to a counterparty is positively correlated with the credit quality of that counterparty. Wrong/right way risk, as an additional source of risk, is rightly of concern to banks and regulators.

Some financial markets are closely interlinked, while others are not. For example, CDS price movements have a feedback effect on the equity market, as a trading strategy commonly employed by banks and other market participants consists of selling a CDS on a reference entity and hedging the resulting credit exposure by shorting the stock. On the other hand, Moody's

Investor 's Service (2000) presents statistics that suggest that the correlations between interest rates and CDS spreads are very small.

To capture wrong/right way risk, we need to determine the dependency between counterparties and to correlate the credit spreads or hazard rates with the other market risk factors, e.g. equities, commodities, etc., in the scenario generation.

We use an equity swap as an example. Assume the correlation between the underlying equity price and the credit quality (hazard rate) of party B is  . The impact of the correlation on the CVA is show in Table 5. The results say that the CVA increases when the absolute value of the negative correlation increases.

Table 5. The impact of wrong way risk on the CVA

This table shows that the CVA increases while the negative correlation  increases in the absolute value. We use an equity swap as an example and assume that there is a negative correlation between the equity price and the credit quality of party B.

Conclusion

This article presents a framework for pricing risky contracts and their CVAs. The model relies on the probability distribution of the default jump rather than the default jump itself, because the default jump is normally inaccessible. We find that the valuation of risky assets and their CVAs, in most situations, has a backward recursive nature and requires a backward induction valuation. An intuitive explanation is that two counterparties implicitly sell each other an option to default when entering into an OTC derivative transaction. If we assume that a default may occur at any time, the default options are American style options. If we assume that a default may only happen on the payment dates, the default options are Bermudan style options. Both

Bermudan and American options require backward induction valuations.

Appendix A. Proofs

Proof of Proposition 1: Under the unilateral credit risk assumption, we only consider the default risk when the asset is in the money. Assume that a default may only occur on the payment date. Therefore, the risky value of the asset at t is the discounted expectation of all possible payoffs and is given by

              t t t F F F T T X T X X X T t F E X T t q T T t D E X T t q T T t p T t D E t V T T T ) , ( ) , ( ) ( 1 1 1 ) , ( 1 ) , ( ) ( ) , ( 1 ) , ( ) ( 0
0 0             (A1a)
where

    ) ( 1 ) , ( 1 1 ) , ( ) , ( 0 T T t q T t D T t F T X      (A1b) Proof of Proposition 2: Let 0 T t  . On the first payment day, let ) ( 1 T V
denote the risky value of the asset excluding the current cash flow 1 X . According to Proposition 1, the risky value of the asset at t is given by

    t F ) ( ) , ( ) ( 1 1 1 0 T V X T T F E t V   (A2a) where       ) ( 1 ) , ( 1 1 ) , ( ) , ( 0 ) ( 1 0 1 0 1 1 T T t q T T D T T F X T V       (A2b) Similarly, we have     1 ) ( ) , ( ) ( T V X T T F E T V F   (A3) Note that ) , ( 1 0 T T F is 1 T F -
                      t t t t t F F F F F F F ) ( ) , ( ) , ( ) ( ) , ( ) , ( ) , ( ) , ( ) ( ) , ( ) ( 2 1 0 1 2 1 1 0 1 2 2 1 2 2 1 1 0 1 1 0 1 1 1 0 1 1 T V T T F E X T T F E T V T T F E X T T F E T T F E X T T F E T V X T T F E t V j j j i i i j j j T T                 (A4)
By recursively deriving from 2

T forward over m T , where

m m X T V  ) ( , we have            m i t i i j j j X T T F E t V 1 1 0 1 ) , ( ) ( F (A5)
Proof of Proposition 3: We assume that a default may only occur on the payment date.

At time T, there are four possible states: 1) both A and B survive, 2) A defaults but B survives, 3)

A survives but B defaults, and 4) both A and B default. The joint distributions of A and B are given by ( 15). Depending on whether the payoff is in the money or out of the money at T, we have

          t T A X B X t T t T AB A A X t T AB B B X X T t k T t k T t D E X T t K E X T t p T T t p T T t p T T t p X T t p T T t p T T t p T T t p T t D E t V T T T T F F F F ) , ( 1 )
                     (A6a) where   ) ( ) ( ) ( 1 ) , ( ) , ( ) , ( ) ( ) , ( ) , ( ) ( ) , ( ) , ( ) ( ) , ( ) , ( ) , ( T T T T t T t q T t q T T t q T t p T T t p T t q T T t p T t p T t k AB B B AB A B AB A B B A B B A B B                (A6b)   ) ( ) ( ) ( 1 ) , ( ) , ( ) , ( ) ( ) , ( ) , ( ) ( ) , ( ) , ( ) ( ) , ( ) , ( ) , ( T T T T t T t q T t q T T t q T t p T T t p T t q T T t p T t p T t k AB A A AB A B AB B A A B A A A B A                (A6c)

B. A practical framework for calculating bilateral CVA

We develop a practical framework for calculating bilateral CVA at counterparty portfolio level based on the theory described above. The framework incorporates netting and margin agreements, and captures right/wrong way risk.

Two parties are denoted as A and B. All calculations are from the perspective of party A.

Let the valuation date be t. The CVA computation procedure consists of the following steps.

B.1. Risk-neutral Monte Carlo scenario generation

One core element of the trading credit risk modeling is the Monte Carlo scenario generation (market evolution). This must be able to run a large number of scenarios for each risk factor with flexibility over parameterization of processes and treatment of correlation between underlying factors. Credit exposure may be calculated under real probability measure, while CVA or pricing counterparty credit risk should be conducted under risk-neutral probability measure.

Due to the extensive computational intensity of pricing counterparty risk, there will inevitably be some compromise of limiting the number of market scenarios (paths) and the number of simulation dates (also called "time buckets" or "time nodes"). The time buckets are normally designed fine-granularity at the short end and coarse-granularity at the far end. The details of scenario generation are beyond the scope of this paper.

B.2. Cash flow generation

For ease of illustration, we choose a vanilla interest rate swap, as interest rate swaps collectively account for around two-thirds of both the notional and market value of all outstanding derivatives (FinPricing (2018)) Cash flow generation for products without early-exercise provision is quite straightforward. For early-exercise products, one can use the approach proposed by [START_REF] Longstaff | Valuing American options by simulation: a simple least-squares approach[END_REF] to obtain the optimal exercise boundaries and then the payoffs.

B3. Aggregation and netting agreements

After generating cash flows for each deal, we need to aggregate them at counterparty portfolio level at each scenario and each time bucket. The cash flows are aggregated by either netting or nonnetting based on the netting agreements. A netting agreement is a provision that allows the offset of settlement payments and receipts on all contracts between two counterparties.

Another important use of netting is the close-out netting that allows the offset of close-out values.

For netting, we add all cash flows together at the same scenario and the same time bucket to recognize offsetting. The aggregated cash flow under netting at scenario j and time bucket k is given by

  i i k j k j , , ,   (B3)
For nonnetting, we divided cash flows into positive and negative groups and add them separately. In other words, the offsetting is not recognized. The aggregated cash flows under nonnetting at scenario j and time bucket k are given by

B4. Margin (or collateral) agreements

For a more detailed discussion on pricing collateralized contract/portfolio, see Xiao (2013b).

B5. CVA Calculation

After aggregating all cash flows via netting, one can price a portfolio in the same manner as pricing a single deal. We assume that the reader is familiar with the least square Monte Carlo valuation model proposed by [START_REF] Longstaff | Valuing American options by simulation: a simple least-squares approach[END_REF] and thus do not repeat some wellknown procedures for brevity.

If the counterparty portfolio is collateralized, we can calculate the risky value based on equation ( 21) of [START_REF] Xiao | An economic examination of collateralization in different financial market[END_REF]. If there is no collateral agreement, we can price the portfolio according to Proposition 4 in this paper.

CVA is by definition the difference between the risk-free portfolio value and the true (or risky or defaultable) portfolio value.

  ) Proof: See the appendix. The risky valuation in Proposition 2 has a backward nature. The intermediate values are vital to determine the final price. For a discrete time interval, the current risky value has a dependence on the future risky value. Only on the final payment date m T , the value of the contract and the maximum amount of information needed to determine the risk-adjusted discount factor are revealed. The coupled valuation behavior allows us to capture wrong/right way risk

  default distributions can be determined by the reduced-form models. The joint distributions of a bivariate Bernoulli variable can be easily obtained via the marginal distributions by introducing extra correlations. Consider a pair of random variables ( A Y , B Y ) that has a bivariate Bernoulli distribution.

  i = 1,…,m. Each cash flow may be positive or negative. The bilateral risky value of the multiple-payment contract is given by Proposition 4: The bilateral risky value of the multiple-payment contract is given by

  where

  -free discount factor based on the interest rate curve simulated at k T .

  Cox arrival process (also known as a doubly stochastic Poisson process) whose first jump occurs at default and is defined as,

					inf	 t	:		t 0	h	(	s	,		s ds )		 	(1)
	where	) h or (t	) h  denotes the stochastic hazard rate or arrival intensity dependent on an , ( t
			t
	exogenous common state

endogenously, but characterizes it exogenously by a jump process. The stopping (or default) time  of a firm is modeled as a t  , and  is a unit exponential random variable independent of t  .

  Two counterparties are denoted as A and B. Let valuation date be t.

					Consider a financial
	contract that promises to pay a	X		0	from party B to party A at maturity date T, and nothing
		T
	before date T. All calculations in the paper are from the perspective of party A. The risk free value
	of the financial contract is given by

  By definition, an

	t F  t



-measurable random variable is a random variable whose value is known at time t t   . Based on the taking out what is known and tower properties of conditional expectation, we have

Table 1 . Payoffs of a bilaterally defaultable contract

 1 This table displays all possible payoffs at time T. In the case of 0 

	X	T	, there are a total of four

T X

Table 3

 3

	shows that if party B has an infinite collateral threshold		

B H

Table 2 . The impact of collateral threshold B H on the CVA This

 2 table shows that given an infinite A H , the CVA increases while B H increases, where B H denotes the collateral threshold of party B and A H denotes the collateral threshold of party A.

	Collateral Threshold B H	10.1 Mil	15.1 Mil	20.1 Mil	Infinite (  )
	CVA	19,550.91	20,528.65	21,368.44	22,059.30

Table 3 . The impact of collateral threshold

 3 This table shows that given an infinite B H , the CVA decreases while denotes the collateral threshold of party A.

	A H increases, where	H	B

A H on the CVA

A H

Table 4 . The impact of the both collateral thresholds on the CVA

 4 denotes the collateral threshold of party A. This reflects the fact that the collaterals have bilateral impacts on the CVA.

	The CVA may increase or decrease while both collateral thresholds change, where	B H denotes
	the collateral threshold of party B and
	Collateral Threshold	H	B	10.1 Mil	15.1 Mil	20.1 Mil	Infinite (  )
	Collateral Threshold	H	A	10.1 Mil	15.1 Mil	20.1 Mil	Infinite (  )
	CVA			25,752.98	22,448.45	23,288.24		22,059.30

A H

  Based on our theory, we propose a novel cash-flow-based framework (see appendix) for calculating bilateral CVA at the counterparty portfolio level. This framework can easily incorporate various credit mitigation techniques, such as netting agreements and margin agreements, and can capture wrong/right way risk. Numerical results show that these credit mitigation techniques and wrong/right way risk have significant impacts on CVA.

  measurable. According to the taking out what is known and tower properties of conditional expectation, we have

An interest rate swaplet

  This figure illustrates the time relationship between an interest rate swaplet and the simulation time buckets. The floating leg of the swaplet is reset at the fixing date f Brownian Bridge or linear interpolation. The linear interpolation is the expectation of the Brownian Bridge. Then we can calculate the payoff of swaplet i at scenario j as , ) and R denotes the fixed rate.The cash flow amount calculated by (B1) is paid on the payment date p i t , . This value should be allocated into the nearest previous time bucket k T as:

	Interest rate curve simulated at	i t ,	f
	Rates
	M T For swaplet i, there are four important dates: the fixing date f T T ,..., , 1 0 ) in each scenario and N cash flows in the sample swap. Let consider scenario j first. i t , , the starting date s i t , , the ending date e i t , and the payment date p i t , . In general, these dates are not coincidently at the simulation time buckets. The time relationship between swaplet i and the simulation time buckets is illustrated in Figure B1. Figure B1: i t , with the starting date s i t , , the ending date e i t , , and the payment date p i t , . The simulation time buckets are 1 1 ,..., ,   k i i T T T . The simulated interest rate curve is starting at f i t , . Both fixed rate payments and floating-rate payments occur on the same payment dates. The cash flow of swaplet i is determined at the fixing date f i t , that is assumed to be between the simulation time buckets j T and 1  j T . First, we need to create an interest rate curve observed at f i t , by interpolating the interest rate curves simulated at j T and 1  j T via either   ) , ( ) , ; ( , , , , , , e i s i e i s i f i i j t t R t t t F N     (B1) where N denotes the notional; ) , ; ( , , , e i s i f i t t t F denotes the simply compounded forward rate reset at f i t , for the forward period ( s i t , , e i t , ); ) , ( , , e i s i t t  denotes the accrual factor or day count fraction for the period ( s i e i t , , T 1  j T k T 1  k T f i t , s i t , e i t , p i t , t j Terms
			22

Assume that party A pays a fixed rate, while party B pays a floating-rate. Assume that there are M time buckets (

Here we use the recovery of market value (RMV) assumption.

2 1 1 T


Marwen Bouasker

Frédéric Grondin

Pierre Mounanga

Arnaud Pertué

Abdelhafid Khelidj

IMPROVED MEASUREMENT METHODS FOR AUTOGENOUS SHRINKAGE OF CEMENT MORTARS AT VERY EARLY AGE

teaching and research institutions in France or abroad, or from public or private research centers.

Improved measurement methods for autogenous shrinkage of cement mortars at very early age

INTRODUCTION

Autogenous shrinkage is known as one of the phenomena involved in premature (micro) cracking of cement-based materials used in Civil Engineering buildings and constructions. In restrained conditions, early age shrinkage can lead to the formation of both micro-cracks, in the surrounding environment of the aggregates, and crossing cracks, in the cementitious matrix. It results in an increase of the material porosity and consequently in a loss of the containment characteristics and an augmentation of the permeability. The durability of concrete buildings depends then, on the dimensional variations of the cementitious matrices at very early age.

Autogenous shrinkage is a physico-chemical phenomenon coupled with mechanical capillary effects due to pore water consumption by the hydration reactions. At 24 hours of hydration, it can represent more than 40% of the autogenous shrinkage value measured after one year of hydration [START_REF] Baroghel-Bouny | Effect of mix-parameters on autogenous deformations of cement pastes -Microstructural interpretations[END_REF].

During the first stages of hydration (since the first contact between water and cement until the setting), the capillary effects can be considered as negligible and autogenous shrinkage is mainly considered as chemical one (also called «Le Chatelier's contraction») [START_REF] Justnes | Correlating the deviation point between external and total chemical shrinkage with the setting time and other characteristics of hydrating cement paste' in 'Shrinkage 2000[END_REF]. Then, the development of the solid mineral skeleton slows down the deformations of the material and the autogenous shrinkage curve starts to diverge from the Le Chatelier's contraction curve. The self-desiccation of the cementitious matrix becomes the driving mechanism of the autogenous shrinkage. The mechanical resistance of the material gradually increases with the progress of the hydration reactions and the interaction between cement paste and aggregates comes into place.

Measurement of cement-based material deformations at very early age is a major problem [START_REF] Sant | Measurement of volume change in cementitious materials at early ages: review of testing protocols and interpretation of results[END_REF]. Indeed, standard methods for shrinkage measurement are suitable to hardened cementitious systems but cannot be applied during the very early age of the material and particularly during the fluid -solid transition (setting period). It is therefore necessary to propose improved experimental methods adapted to the complex behaviour of fresh cementbased matrices.

The present research work forms part of an extensive study on the influence of matrixinclusion interactions on autogenous shrinkage of mortars. Firstly, five measurement methods especially adapted for the study of the early age deformation of cementitious matrices are described. These methods are then used to measure the volumetric and linear variations of hydrating cement pastes and mortars. The effects of granular inclusions type, water-to-cement ratio (W/C) and curing temperature are analyzed and discussed.

EXPERIMENTAL PROGRAM

Five different measurement methods have been developed for the study of the early age volume changes of cement-based materials: these experimental methods permit to measure continuously and automatically, from the end of casting to 48 hours, the chemical shrinkage, the volumetric autogenous shrinkage, the linear autogenous shrinkage in both horizontal and vertical directions and the restrained autogenous deformations leading to the premature cracking of the material.

These methods are applied to cement pastes with different W/C ratios (W/C= 0.25; 0.30 and 0.40) and mortars prepared with various types of granular inclusions (glass balls and natural sand). The cement used is a Type I ordinary Portland cement (CEM1).

During the tests, precautions are taken to maintain quasi-isothermal conditions in order to decouple pure physico-chemical deformations from thermal deformations caused by the hydration heat release. The experimental systems used for free volume or length change recording are placed in a thermostatic water bath during the whole period measurement (see figure 1). The ring test device used for the measurement of restrained deformations is equipped with a peripheral water circulation system acting as a thermal regulation.

The following sections describe more precisely each experimental apparatus.

Chemical shrinkage measurement

The chemical shrinkage is measured ensuring the permanent water saturation of the hydrating material pore space. A gravimetric method is used [START_REF] Boivin | Experimental assessment of chemical shrinkage of hydrating cement paste[END_REF][START_REF] Holt | Early age autogenous shrinkage of concrete[END_REF]: it consists in recording the buoyancy variation of a sample of cement paste or mortar placed in a flask with an orifice on its cap. The whole system is immersed in water and hung to a balance. The detailed description of the method and protocol can be found in [START_REF] Mounanga | Predicting Ca(OH) 2 content and chemical shrinkage of hydrating cement pastes using analytical approach[END_REF]. The main parameters of the test method are given in table 1.

2.2

Dynamic volumetric autogenous shrinkage measurement This method is based on the "condom method", which consists in recording the buoyancy variation of a latex membrane filled with fresh material, immersed in a liquid bath and hung to a balance. Several researchers already pointed out the experimental artefacts of this test method [START_REF] Justnes | Total and external chemical shrinkage of low W/C ratio cement pastes[END_REF][START_REF] Barcelo | Linear vs. volumetric autogenous shrinkage measurement: material behaviour or experimental artefact?[END_REF]. Among them, the pressure of the membrane on the solidifying matrix, the bleeding phenomenon and the water absorption by the membrane are the most significant.

In order to limit water absorption by the latex membrane, the empty membrane is saturated during 24 hours before starting the test. The measurement system has also been improved to eliminate the effect of bleeding by putting the sample in rotation during the whole test period. The device, presented in figure 2, is composed of a PVC cylinder equipped with a toothed wheel connected to a micro-motor by the intermediary of a driving belt. The micro-motor has been waterproofed with silicon gel. Its rotation speed is controlled by a regulated power supply. The micro-motor, the PVC cylinder and the transmission system (toothed wheel and belt) are assembled on a PVC base hung, by nylon threads, to a balance linked to a computer for automatic recording of the system buoyancy variations. The measurement of volume changes starts ten minutes after the placing of the sample in the PVC cylinder. Table 2 recapitulates the principal parameters of the testing method.

Horizontal and vertical linear autogenous shrinkage measurement

The linear autogenous deformations of cementitious matrices are measured in both horizontal and vertical directions (see figure 3). The experimental system is composed of a ringed flexible moulded polyvinyl chloride (PVC) membrane filled with fresh cement paste or mortar. Precautions are taken to avoid the introduction of air bubbles during the casting of the material. The membrane is sealed with metallic caps and positioned on a rigid PVC frame.

For the "horizontal device", two sensors are placed at each end of the frame. The sensors are waterproof non-contacting displacement measuring system operating on the eddy current principle. This characteristic (waterproof sensors) is important because it makes it possible to maintain the experimental device under water and keep the sample in quasi-isothermal conditions.

For the "vertical test" device, the length changes are measured by a laser sensor located at the top of the experimental system. This sensor is not waterproof: only the lower part of the frame and the membrane filled with cement-based material are immersed in water during the test.

Restrained deformations: quasi-isothermal ring test device

The deformations in restrained conditions are measured using the "ring test" method [START_REF] Hossain | The role of specimen geometry and boundary conditions on stress development and cracking in the restrained ring test[END_REF]. The quasi-isothermal conditions during the tests are ensured by a peripheral water circulation system, whose temperature is controlled by a thermostatic water bath. Photographs and dimensions of the ring test device are given on figures 4 and 5. At the beginning of each test, the fresh cement-based material is cast around the steel ring equipped with strain gauges connected to a data acquisition system.

The device permits to follow automatically the deformations of the steel ring and to detect the moment of the cement-based material cracking. In this study, two isothermal curing temperatures (T=20°C and 25°C) have been imposed to the samples.

RESULTS

3.1

Chemical shrinkage Figure 6 presents the evolution of the chemical shrinkage between 1 and 24 h of hydration of mortars prepared with different W/C ratios (0.30 and 0.40), aggregates of different types (natural sand and glass balls), diameters (1 mm and 2 mm) and concentrations A/C (A/C=0.5 and 1). The curves obtained are almost superimposed up to 8 h of hydration; beyond this first period, they form a tight spindle until the end of the investigation period.

A modification of W/C ratio, granular concentration or type only causes a weak variation of the kinetics and the amplitude of the chemical shrinkage of the studied mortars. A 24 h, the maximum difference observed, equal to 11%, is located between the final curve amplitudes of the mortars with the same W/C and A/C prepared with aggregates of different sizes ("Sand (1mm) W/C=0.4 A/C=0.5" and "Sand (2mm) W/C=0.4 A/C=0.5"). This difference between the two curves can be explained by a more efficient mixing in the presence of the "1 mm" aggregates. Indeed, considering a constant granular mass fraction, the finer aggregates are more numerous and their granular specific surface per gram of sample is larger: the number of cement clusters broken by shearing or impact with sand grains is then larger and the amount of cement directly in contact with water becomes more important. It results in higher amplitude of chemical shrinkage of the mortar prepared with "1 mm" aggregates.

3.2

Volumetric autogenous deformation It is important to remind the complexity of the measurement of cementitious matrix volume changes in particular for the matrices with high W/C ratio because of the bleeding risk. Justnes et al. [START_REF] Justnes | Total and external chemical shrinkage of low W/C ratio cement pastes[END_REF] already observed that this phenomenon tends to over-estimate the autogenous deformations at very early age.

In order to avoid this artefact, a rotary volume change measurement device has been developed: it enables to measure the volume variations of a latex membrane filled with cement paste or mortar and put in rotation during the whole investigation period. For the validation of the system, we applied it to two different cement pastes, the first without bleeding (W/C=0.25) and the second presenting a strong bleeding risk (W/C=0.4).

Figure 8 displays the results obtained with the W/C=0.25-cement paste. The curves are initialized at 1 h of hydration in order to eliminate the disturbances observed at the beginning of the test. These disturbances are due to the movement of the sample, which is still fluid, inside the cylinder in rotation. The comparison between the curves obtained with the rotary device and those measured by Mounanga [START_REF] Mounanga | Experimental study on the behaviour of cement pastes at very early age: hydration, shrinkage and thermal properties[END_REF] in static conditions (no rotation) shows that the both test methods give very similar results for cement pastes with low W/C ratio.

Figure 9 presents the results of volumetric autogenous shrinkage measurement of cement pastes with W/C=0.40 with and without rotation of the sample. These results are compared with the chemical shrinkage evolution of W/C=0.40-cement pastes presented by Mounanga et al. in 2004 [6]. It can be noted that in static conditions, the period of volume changes corresponding to Le Chatelier's contraction is much longer than in the case of dynamic measurement. This phenomenon is explained by the formation of a water film between the membrane and the sample in static conditions which artificially extends the period of Le Chatelier's contraction. During the dynamic test, this film does not have time to be formed and the bleeding effect is avoided. The initializing of the deformations at 15 hours of hydration (figure 10), i.e. beyond the knee-point marking the transition between chemical deformations and self-desiccation shrinkage, underlines that the volume changes measured in static and in dynamic conditions are completely superimposed: at this stage of hydration, the bleeding of the sample does not have any more influence on the material deformation.

As a matter of fact, for cement-based systems with high W/C ratio, measurement methods in static conditions or with continuous rotation cannot provide the real autogenous deformations observed in actual field conditions:

-In static conditions, the measured shrinkage is very over-estimated due to the extension of Le Chatelier's contraction phase induced by the bleeding phenomenon; -In continuous rotation, bleeding is avoided and the sample remains homogeneous during the test period. Field conditions are situated somewhere between these two test conditions. Therefore, the function of the dynamic volumetric autogenous shrinkage device is not to reproduce completely field conditions but rather to permit the measurement of the volume change of a "perfectly" homogeneous sample, decoupling physico-chemical phenomena (chemical shrinkage, setting, self-desiccation) from bleeding and sedimentation phenomena.

Uniaxial autogenous deformation

The measurement of the volume deformations (chemical or volumetric autogenous shrinkage) does not provide sufficient information to describe the material behaviour. Indeed, for concrete structure modelling, it is important to know the uniaxial deformations of the material. An improved horizontal linear autogenous shrinkage device has been developed and a second device for vertical measurement is under development.

The first results obtained with the horizontal device are presented in figures 11 and 12. One can notice in figure 11 a very good repeatability of the signal given by the two sensors located on both ends of the sample. The sensors are waterproof non-contact eddy current transducers. It is thus possible to carry out the tests in thermostatic water bath in order to maintain the sample at constant temperature. After 20 minutes of test, the sample temperature is stabilized at 20±1°C. The test method is used to investigate the effect of W/C ratio on the evolution of autogenous shrinkage. Figure 12 shows the lower the W/C ratio the higher the autogenous shrinkage. This result is in agreement with previous studies on the subject [START_REF] Baroghel-Bouny | Effect of mix-parameters on autogenous deformations of cement pastes -Microstructural interpretations[END_REF]: in a cement paste with low water content, the self-desiccation is faster and more intense and thus causes a more important shrinkage. The periodic disturbances observed on figure 12

3.4

Isothermal ring test This ring test device was developed to quantify, in autogenous, quasi-isothermal and restrained conditions, the deformations of cement-based systems as well as the age of cracking, from the measuring of the steel ring deformations. Figure 13 displays the deformations of the steel ring until the cracking of the cement paste (W/C=0.3) at 20°C and 25°C. The result analysis shows that the temperature increasing causes an augmentation of the deformation rate and a shortening of the age of cracking. This phenomenon is explained by the thermoactivated characteristic of hydration reactions, which induces an acceleration of both deformations and internal stresses when the temperature increases.

In order to underline the kinetic effect of temperature, it is preferable to plot the deformations as a function of hydration degree. A continuous function of the hydration degree has been obtained using the NIST's model CEMHYD3D [START_REF] Bentz | CEMHYD3D: A three-dimensional cement hydration and microstructure development modelling package[END_REF]: the value of the β-parameter, which enables to convert the number of model cycles to real time, has been fitted from experimental values measured at different temperatures and water-to-cement ratios (figure 14). In our case, the value β=0.00014 gives the best fitting.

Numerical results (CEMHYD3D)

W/C=0.3 T=20°C W/C=0.4 T=20°C W/C=0.3 T=10°C (ref. [START_REF] Mounanga | Predicting Ca(OH) 2 content and chemical shrinkage of hydrating cement pastes using analytical approach[END_REF]) W/C=0.3 T=20°C (ref. [START_REF] Mounanga | Predicting Ca(OH) 2 content and chemical shrinkage of hydrating cement pastes using analytical approach[END_REF]) W/C=0.3 T=30°C (ref. [START_REF] Mounanga | Predicting Ca(OH) 2 content and chemical shrinkage of hydrating cement pastes using analytical approach[END_REF]) W/C=0.4 T=10°C (ref. [START_REF] Mounanga | Predicting Ca(OH) 2 content and chemical shrinkage of hydrating cement pastes using analytical approach[END_REF]) W/C=0.4 T=20°C (ref. [START_REF] Mounanga | Predicting Ca(OH) 2 content and chemical shrinkage of hydrating cement pastes using analytical approach[END_REF]) W/C=0.4 T=30°C (ref. [START_REF] Mounanga | Predicting Ca(OH) 2 content and chemical shrinkage of hydrating cement pastes using analytical approach[END_REF]) Coefficient β=0.00014

CONCLUSIONS

During this study, various experimental techniques for the measurement of free and restrained deformations of cementitious matrices are described and the first results obtained are analyzed. It was in particular shown that:

The presence of granular inclusions slightly accelerates the evolution of mortar chemical shrinkage: the crossing analysis of the results highlighted that the additional destruction of cement grain clusters during mixing can be considered as the main cause of Le Chatelier's contraction acceleration at very early age. This phenomenon does not influence significantly the quasi-linear relation between chemical shrinkage and hydration degree of the mortars. The bleeding phenomenon influences the measurement of cement paste volume change at very early age, by inducing an extension of the deformation phase of autogenous shrinkage dominated by Le Chatelier's contraction. This effect, which has be found significant for cement pastes with W/C ratio≥0.4, leads to an overestimating of autogenous shrinkage. The rotary test device developed in the present study makes it possible to eliminate this artefact.

The use of non-contact and waterproof eddy current sensors enables to quantify precisely, in quasi-isothermal conditions and since the casting of the material, the uniaxial autogenous deformations of the cementitious systems. The age of cracking measured with the ring test device decreases when the temperature increases. This result can be explained by the acceleration of hydration reactions and of the development of internal stresses. However, in the temperature range investigated (20-25°C), cracking appears at nearly the same cement hydration degree (≈58±2% of hydration).

Figure 2 :

 2 Figure 2: Dynamic volumetric autogenous shrinkage measurement device (a-Sample, b-Rotary system)

Figure 3 :

 3 Figure 3: Linear autogenous shrinkage measurement devices (a-Horizontal system, b-Waterproof eddy-current sensors, c-Flexible PVC membrane, d-Vertical system)

Figure 4 :Figure 5 :

 45 Figure 4: Diagrams (a-Top view, b-Side view) and general view (c) of the quasi-isothermal ring test device

5 Figure 6 :Figure 7 :

 567 Figure 6: Chemical shrinkage vs. hydration age at T=20 ± 1°C

Figure 8 :

 8 Figure 8: Volumetric autogenous shrinkage vs. hydration age (Cement paste with W/C=0.25 at T=20 ± 1°C)

Figure 9 :Figure 10 :

 910 Figure 9: Comparison between chemical shrinkage and volumetric autogenous shrinkage measured with and without rotation of the sample (Cement paste with W/C=0.4 at T=20 ± 1°C)

Figure 11 :

 11 Figure 11: Horizontal linear autogenous shrinkage vs. hydration age (Cement paste with W/C=0.3 at T=20 ± 1°C)

3 Figure 12 :

 312 Figure 12: Effect of the W/C ratio on the horizontal linear autogenous shrinkage (Cement pastes with W/C=0.25 and 0.3 at T=20 ± 1°C)

Figure 13 :

 13 Figure 13: Steel ring deformations vs. hydration age (Cement paste with W/C=0.3)

Figure 14 :Figure 15 :

 1415 Figure 14: Comparison between numerical and experimental results of hydration degree of CEM1 cement pastes at different temperatures

Table 1 :

 1 Test parameters for Le Chatelier's contraction measurement

	Test parameters	Values
	Sample mass	20 -30 g
	Exchange surface area between sample and water	20 ± 1 cm 2
	Sample thickness	9 ± 1 mm
	Liquid of immersion	Distilled water
	Water bath temperature	20 ± 1°C

Table 2 :

 2 Test parameters for dynamic volumetric autogenous shrinkage measurement

ACKNOWLEDGEMENTS

The authors gratefully acknowledge the support received from the Chamber of Commerce of Nantes-Saint-Nazaire and Nantes Metropole, which actively participate in funding the research work of the first author. Thanks to Ahmed Loukili and Philippe Turcry for fruitful discussion.


Ashraf Ahmed
email: ashraf.ahmed@univ-nantes.fr

Abdul-Hamid Soubra
email: abed.soubra@univ-nantes.fr

Improved Subset Simulation for the SLS Analysis of Two Neighboring Strip Footings Resting on a Spatially Random Soil

Keywords: subset simulation, improved subset simulation, conditional simulation, strip footings, differential settlement

 to efficiently calculate the small failure probability. In the present paper, a more efficient approach called the improved Subset Simulation (iSS) is employed. In this approach the efficiency of SS is increased by replacing the first step of SS by a conditional simulation in which the realizations are generated outside a hypersphere of a given radius. This approach is illustrated here through the probabilistic analysis at the serviceability limit state (SLS) of two neighboring strip footings that rest on a soil with spatially varying Young's modulus. A comparison between SS and iSS approaches has shown that a considerable reduction in the number of realizations can be achieved when using the iSS approach.

INTRODUCTION

The classical Monte Carlo Simulation (MCS) methodology is generally used to calculate the failure probability of problems involving a spatial variability of the soil properties. This method is very time-consuming when computing a small failure probability. This is due to the large number of realizations required in such a case. As alternative to MCS methodology, the Subset Simulation (SS) approach was proposed by [START_REF] Au | Estimation of small failure probabilities in high dimentions by subset simulation[END_REF] to calculate the small failure probability. The first step of SS method is to generate a given number of realizations of the uncertain parameters using the classical MCS technique. The second step is to use the Metropolis-Hastings (M-H) algorithm to generate realizations in the direction of the limit state surface (i.e. G=0). This step is repeated until reaching the limit state surface. It should be emphasized here that in case of a small failure probability, SS requires the repetition of the second step many times to reach the limit state surface. This leads to a high computational time. To reduce the computational cost of SS, [START_REF] Defaux | Mécanique probabiliste et intégrité des cuves des réacteurs nucléaires[END_REF] proposed an improved subset simulation (iSS) method. In this method, the efficiency of SS is increased by replacing the first step by a conditional simulation. In other words, instead of generating realizations directly around the origin by the classical MCS, the realizations are generated outside a hypersphere of a given radius. Consequently, the number of realizations required to reach the limit state surface is significantly reduced. Notice that [START_REF] Defaux | Mécanique probabiliste et intégrité des cuves des réacteurs nucléaires[END_REF] have employed the iSS to calculate the failure probability in the case where the uncertain parameters are modeled by random variables. In the present paper, the iSS is employed in the case where the uncertain parameters are modeled by random fields. This method is illustrated through the computation of the probability (P e ) of exceeding a tolerable differential settlement between two neighboring strip footings resting on a soil with a spatially varying Young's modulus. The footings are subjected to central vertical loads with equal magnitude. The random field is discretized using the Karhunen-Loeve (K-L) expansion. The differential settlement between the two footings was used to represent the system response. The deterministic model used to compute the system response is based on numerical simulations using the commercial software FLAC.

REVIEW OF THE CLASSICAL SUBSET SIMULATION (SS) APPROACH

Subset simulation was proposed by [START_REF] Au | Estimation of small failure probabilities in high dimentions by subset simulation[END_REF] to compute the small failure probabilities. The basic idea of the SS approach is that the small failure probability can be expressed as a product of larger conditional failure probabilities. Consider a failure region F defined by the condition G<0 where G is the performance function and let (s 1 , …, s k , ..., s Nt ) be N t samples located in the space of the uncertain variables where s represents a vector of random variables. It is possible to define a sequence of nested failure regions F 1 , …, F j , ..., F m of decreasing size where F

F ... F ... F m j 1 = ⊃ ⊃ ⊃ ⊃
as shown in figure 1. An intermediate failure region F j can be defined by G j <C j where C j >0. Thus, there is a decreasing sequence of positive numbers C 1 , …, C j , ..., C m corresponding respectively to F 1 , …, F j ,…, F m where C 1 >…>C j >...> C m =0. The N t samples (s 1 , …, s k , ..., s Nt ) will be divided into groups of equal number N s of samples (s 1 , …, s k , ..., s Ns ). Thus, N t =mxN s where m is the number of failure regions. The Ν s samples of the first group are generated by MCS methodology according to a target PDF (P t ). The corresponding failure probability P(F 1 ) is calculated as follows:

∑ = = Ns 1 k k F s 1 ) s ( I N 1 ) F ( P 1 (1) where 1 I 1 F = if 1 F s ∈ and 0 I 1 F = otherwise.
On the other hand, the samples of the remaining groups are generated using Metropolis-Hastings algorithm according to a proposal PDF (P p ). The failure probability corresponding to an intermediate failure region F j where j≠1 is calculated as follows:

∑ = - = Ns 1 k k F s 1 j j ) s ( I N 1 ) F F ( P j (2) 1 I j F = if

IMPROVED SUBSET SIMULATION (iSS) APPROACH

As mentioned previously, the basic idea of iSS is to replace the first step of SS by a conditional simulation [START_REF] Defaux | Mécanique probabiliste et intégrité des cuves des réacteurs nucléaires[END_REF]) in which, the realizations are generated outside a hypersphere of a given radius β (figure 2). Based on this conditional simulation, the failure probability P(F 1 ) corresponding to the first level of subset simulation (i.e. level 0) is calculated as follows [START_REF] Defaux | Mécanique probabiliste et intégrité des cuves des réacteurs nucléaires[END_REF]:

∑ = β χ - = Ns 1 k k F s 2 M 1 ) s ( I N 1 )) ( 1 ( ) F ( P 1 (4)
where M χ is the chi-square distribution with M degrees of freedom and 1 I 1 F = if 1 F s ∈ and 0 I 1 F = otherwise. The advantage of using the conditional simulation is to generate realizations in the proximity of the limit state surface leading to a reduction in the number of realizations required to reach this surface. Notice finally that the realizations of the remaining levels (i.e. levels 1 to m-1) are generated using Metropolis-Hastings algorithm. The failure probability of a given level j where j≠1 is calculated using Eq. ( 2) and the final failure probability P(F) is calculated using Eq.

(3).

IMPLEMENTATION OF iSS IN THE CASE OF RANDOM FIELDS

As mentioned before, this paper aims at employing the iSS approach for the computation of the failure probability in the case of a spatially varying soil property. To achieve this purpose, a link between SS and K-L expansion is performed. For a Gaussian random field E(X, θ), where X denotes the spatial coordinates and θ indicates the random nature of this random field, the random field can be approximated by the K-L expansion as follows [START_REF] Spanos | Stochastic Finite Element Expansion for Random Media[END_REF]:

) ( ) X ( ) (X, E i i M 1 i i θ ξ φ λ + μ ≈ θ ∑ = (5)
where μ is the mean of the random field, M is the size of the series expansion, λ i and i φ are the eigenvalues and eigenfunctions of the covariance function C(X 1 , X 2 ), and ξ i (θ) is a vector of standard uncorrelated random variables. In the present paper, the random field E was assumed to follow a log-normal probability density function so that ln(E) is a normal random field with mean value μ ln and standard deviation σ ln . For a lognormal random field, the K-L expansion given in Eq. ( 5) becomes [START_REF] Cho | Effect of spatial variability of cross-correlated soil properties on bearing capacity of strip footing[END_REF]:

⎥ ⎦ ⎤ ⎢ ⎣ ⎡ θ ξ φ λ + μ ≈ θ ∑ = ) ( ) X ( exp ) (X, E i i M 1 i i ln (6)
On the other hand, the random field was assumed to follow an exponential covariance function as follows:

⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ - - - - σ = y ln 2 1 x ln 2 1 2 ln 2 2 1 1 l y y l x x exp )] y , x ( ), y , x [( C (7)
where (x 1 , y 1 ) and (x 2 , y 2 ) are the coordinates of two arbitrary points in the domain over which the random field is defined and l ln x and l ln y are respectively the horizontal and vertical lengths over which the values of ln(E) are highly correlated.

Notice that in the case of an exponential covariance function, the eigenvalues and eigenfunctions are given analytically. Their solutions are presented in [START_REF] Spanos | Stochastic Finite Element Expansion for Random Media[END_REF].

Notice that the basic idea of the link between SS and K-L expansion was given in [START_REF] Ahmed | Subset simulation and its application to a spatially random soil[END_REF] and is briefly presented herein in the case of iSS approach applied to a random field problem. This link is performed through the standard normal random variables

M ..., , 1 i i } { = ξ
appearing in Eq. ( 6) as follows: for a given random field realization discretized by K-L expansion, the system response is calculated in two steps. The fist step is to substitute the vector

M ..., , 1 i i } { = ξ
of this realization in Eq. ( 6) to calculate the value of the random field at each point in the domain according to its coordinates. The second step is to use the deterministic model to calculate the corresponding system response. The algorithm of iSS proposed in this paper for the case of a spatially varying soil property can thus be described as follows:

1. Prescribe a radius β for the hypersphere and generate a vector of standard normal random variables {ξ 1 , …, ξ i , ..., ξ M } by MCS methodology. This vector must realize that its norm is larger than the prescribed radius β. 2. Substitute the vector {ξ 1 , …, ξ i , ..., ξ M } in the K-L expansion (Eq. 6) to obtain the first realization of the random field. Then, use the deterministic model to calculate the corresponding system response. 3. Repeat steps 1 and 2 until obtaining a prescribed number N s of realizations and their corresponding system response values. Then, evaluate the corresponding values of the performance function to obtain the vector

} G ,..., G ,..., G { G Ns 0 k 0 1 0 0 =
. Notice that the subscript '0' refers to the first level (level 0). 4. Prescribe a constant intermediate failure probability P(F j ) for all the failure regions F j and evaluate the first failure threshold C 1 which corresponds to the failure region F 1 where C 1 is equal to the [(N s xP(F j ))+1] th value in the increasing list of elements of the vector G 0 . Thus, among the N s realizations, there are [N s xP(F j )] ones whose values of the performance function are less than C 1 (i.e. they are located in the region F 1 ). The failure probability P(F 1 ) is then calculated by Eq. ( 4) 5. The vectors of {ξ 1 , …, ξ i , ..., ξ M } corresponding to the realizations that are located in the region F 1 (from step 4) are used as 'mother vectors' to generate additional [(1-P(F j ))N s ] vectors of {ξ 1 , …, ξ i , ..., ξ M } using Metropolis-Hastings algorithm. These new vectors are substituted in Eq. ( 6) to obtain the random field realizations of level 1. Then, the values of the performance function corresponding to these realizations should be calculated and gathered in an increasing order in the vector of performance functions

} G ,..., G ,..., G { G Ns 1 k 1 1 1 1 = . 6.

EXAMPLE PROBLEM

In this section, the efficiency of the iSS is illustrated through an example problem. In this example, a probabilistic analysis at SLS of two neighboring strip footings resting on a soil with a spatially varying Young's modulus was performed. Each footing is subjected to a central vertical load P=1000kN/m. Due to the soil heterogeneity, the two footings exhibit a differential settlement δ. This differential settlement is calculated as the absolute difference

2 1 δ - δ = δ
where δ 1 and δ 2 are the settlements (computed at the footing center) of the two footings. The differential settlement δ is used to represent the system response. The Young's modulus is modeled by a random field and it was assumed to follow a log-normal distribution.

Its mean value and coefficient of variation are respectively μ E =60MPa and COV E =15%. It was discretized using K-L expansion. An exponential covariance function (Eq. 7) was used in this paper. Although an isotropic random field is often assumed in literature (e.g. Fenton andGriffiths 2002, 2005), the vertical autocorrelation length tends to be shorter than the horizontal one. A common ratio of about 1 to 10 for these autocorrelation lengths can be used [START_REF] Baecher | Reliability and Statistics in Geotechnical Engineering[END_REF]. Notice however that in this paper, different values of autocorrelation lengths were studied and analyzed in order to explore some interesting features related to the autocorrelation lengths. The performance function used to calculate the probability (P e ) of exceeding a tolerable differential settlement is defined as follows:

δ - δ = max G (9)
where δ max is a prescribed tolerable differential settlement and δ is the computed differential settlement due to the soil heterogeneity. The deterministic model used to calculate the differential settlement δ is based on numerical simulations using FLAC.

For this computation, two footings (each of width b) were considered in the analysis (Figure 3). The two footing centers are separated by a distance D. A non-uniform optimal mesh composed of 1290 zones was employed. Although this paper presents an SLS analysis, the soil behavior was modeled by a conventional elastic-perfectly plastic model based on Mohr-Coulomb failure criterion in order to take into account the possible plastification that may take place near the footing edges even under the service loads. The strip footings were modeled by a linear elastic model. They are connected to the soil via interface elements. The values of the different parameters of the soil, footings and interfaces are given in Table 1. In order to calculate the differential settlement for a given random field realization, (i) the coordinates of the center of each element of the mesh were calculated; then, Eq. ( 6) was used to calculate the value of the Young's modulus at each element, (ii) geostatic stresses were applied to the soil, (iii) the obtained displacements were set to zero in order to obtain the footing displacement due to only the footings applied loads and finally, (iv) the service loads were applied to the footings and the vertical displacements at the footings centers due to these loads are calculated. The differential settlement is calculated as the absolute difference between the two footings displacements. It should be noticed that, for all the probabilistic analyses performed in this paper, a normal PDF was used as a target PDF (i.e. it was used to generate the N s realizations of the first level). However, a uniform PDF was chosen as a proposal PDF (i.e. it was used to generate the N s realizations of each one of the remaining levels). The intermediate failure probability P(F i ) was chosen equal to 0.1. Also, the tolerable differential settlement δ max was assumed equal to 3.5x10 -3 m. It should be mentioned here that l ln x and l ln y were normalized with respect to the distance D between the centers of the two footings (i.e. L ln x =l ln x /D and L ln y =l ln y /D). The numerical results have shown that this assumption is valid when the ratio D/b is constant. Notice that all the probabilistic results presented in this paper correspond to a ratio of D/b=2. Finally, the number of terms of K-L expansion used in this paper is M=100 terms except in the two following cases (L ln x =L lny =0.25 in case of an isotropic random field and, L ln x ≤0.25 and L ln y ≤0.25 in case of an anisotropic random field) where 500 terms were required.

Probabilistic results

This section presents the probabilistic results of the example problem described above. First, the number of realizations N s to be used in the different levels of iSS has to be selected. This number must be sufficient to provide accurate P e values. Different values of N s were considered. For each N s value the failure thresholds C 1 , C 2 , etc. were calculated and presented in Table 2 when the radius β is equal to zero. This table shows that the value of the failure threshold decreases with the successive levels until reaching a negative value at the last level. -0.00009 -0.00051 -0.00036 -0.00037 -0.00039 -0.00038 For each N s value presented in Table 2, P e corresponding to each level j was calculated by iSS as follows:

) ( ... ) ( ) ( ) ( 1 1 2 1 - = j j j F F P x x F F P x F P F P ( 11
)
These P e values were compared to those computed by the crude MCS methodology using a number N=30,000 realizations (Figure 4). Notice that at a given level j, the P e value is calculated by MCS methodology as follows:

∑ = = N k k F j G I N F P 1 ) ( 1 ) ( (12)
in which, G k is the value of the performance function at the k th realization and I F =1 if G k <C j and I F =0 otherwise. The comparison has shown that for N s =1,000 realizations, the P e value computed by iSS at the different levels is very close to that computed by the crude MCS methodology. Thus N s =1,000 realizations will be used in all the probabilistic analyses performed in this paper. Notice that when N s =1,000 realizations, the coefficient of variation of P e by iSS is COV Pe =31.5%. However, COV Pe by MCS methodology using 30,000 realizations is equal to 31.3%. Notice also that when β=0, four levels of SS were required to reach the limit state surface G=0. This means that a total number of realization N t =1,000+(900x3)=3,700 realizations were required to calculate P e . Thus, for the same accuracy, the number of realizations (and consequently, the computation time) is reduced by 87.7% with respect to MCS when β=0. (i.e. when the classical SS is used). This number can be again reduced by increasing β (i.e. by using iSS). Table 3 shows that, when β increases, the total number of realizations decreases. When β=11.5, only 2 levels are required. Thus, the total number of realizations is N t =1000+900=1,900 realizations. As a conclusion, the number of realizations (and consequently, the computation time) required by SS approach could be reduced by 48.6% by employing the iSS approach.

Effect of the autocorrelation length on P e in the case of an isotropic random field

Figure 5 shows the effect of the autocorrelation length on the P e value in the case of an isotropic random field. This figure indicates that P e presents a maximum value when L ln x =L ln y =1. This can be explained by the fact that when the autocorrelation lengths are very small, one obtains a highly heterogeneous soil in both the vertical and the horizontal directions with a great variety of high and small values of the Young's modulus beneath the footings. In this case, the soil under the footings contains a mixture of stiff and soft soil zones. Due to the high rigidity of the footings, their movements are resisted by the stiff soil zones. This leads to a small value of the footings displacements (i.e. to a small differential settlement) and thus, to a small value of P e . On the other hand, when the autocorrelation lengths are large, the soil tends to be homogenous. This means that the differential settlement tends to be very small (close to zero) which leads to a very small value for P e . For the intermediate values of the autocorrelation lengths corresponding to (P e ) max , there is a high probability that one footing rests on a stiff soil zone and the other on a relatively soft soil zone. This leads to a high differential settlement and thus to a high P e value. This configuration corresponds to the case where L ln x =L ln y =1.

Effect of the autocorrelation lengths on P e in the case of anisotropic random field

Figures 6 shows the effect of L ln x on P e when L ln y =0.25. This figure shows that P e presents a maximum value when L ln x =1. For the very small values of L ln x compared to L ln y , one obtains a vertical multilayer composed of thin sub-layers where each sub-layer may have a high or a small value of the Young's modulus. The sub-layers with high values of the Young's modulus prevent the movements of both footings and thus lead to a small value of P e . On the other hand, when L ln x is very large compared to L ln y , one obtains a horizontal multilayer (case of a onedimensional vertical random field) for which each sub-layer may have a high or a small value of the Young's modulus. This leads to the same displacement for both footings and thus to a very small value of P e . Finally, when L ln x is equal to 1, the horizontally extended stiff layers become less extended and thus one obtains a greater differential settlement and consequently a greater value of P e . The effect of L ln y is presented in Figure 7 when L ln x =2.5. This figure shows that the P e value increases with the increase in L ln y . This can be explained as follows: when L ln y is very small, the two footings rest on a horizontal multilayer composed of thin sub-layers where each sub-layer may have a high or a small value of the Young's modulus. This means that δ 1 and δ 2 are almost equal. Thus, the differential settlement δ is very small which results in a small value of P e . On the other hand, when L ln y is very large, the soil tends to the case of a one-dimensional horizontal random field. In this case, one obtains vertically extended stiff sub-layers adjacent to vertically extended soft sub-layers. For the chosen value of L ln x , there is a high probability that one footing rests on a vertical stiff layer and the other one rests on a vertical soft layer which leads to a high differential settlement and thus to a great value of P e .

CONCLUSION

This paper presents an efficient method to perform a probabilistic analysis of geotechnical structures that involve spatial variability. This method is an improvement of the classical subset simulation approach to calculate the small failure probabilities using a reduced number of realizations. It was illustrated through an example problem in which, a probabilistic analysis at SLS of two neighboring strip footings was performed. The footings rest on a soil with a spatially varying Young's modulus. The proposed procedure has significantly reduced the number of realizations required by the classical subset simulation approach to calculate the probability P e of exceeding a tolerable differential settlement. A parametric study to investigate the effect of the autocorrelations lengths on P e has shown that: (i) in case of an isotropic random field, P e presents a maximum value when the autocorrelation lengths are equal to the distance between the footings centers, (ii) in case of an anisotropic random field, for a given L ln y value, P e presents a maximum when L ln x =D. However, for a given value of L ln x , P e increases with the increase in L ln y and then it attains an asymptote corresponding to the case of a one-dimensional random field.

F

  = otherwise. The failure probability P(F) of the failure region F can be calculated from the sequence of conditional failure probabilities as follows:

Figure

  Figure 1. Nested failure domain Figure 2. Samples generation outside a hypersphere of radius β

  Evaluate the second failure threshold C 2 as the [(N s xP(F j ))+1] th value in the increasing list of the vector G 1 . 7. Repeat steps 5 and 6 to evaluate the failure thresholds C 3 , C 4 , …, C m corresponding to the failure regions F 3 , F 4 , …, F m . Notice that contrary to all other thresholds, the last threshold C m is negative. Thus, C m is set to zero and the conditional failure probability of the last level [P(F m |F m -1)] is calculated as: , the failure probability P(F) is evaluated according to Eq. (3).

Figure 3 .

 3 Figure 3. Soil domain and mesh used in the numerical simulations

Figure 4 .

 4 Figure 4. Comparison between P e computed by iSS and that computed by MCS methodology at each level of iSS (β=0, L ln x =2.5 and L ln y =0.25).

Figure 5 .

 5 Figure 5. Effect of the autocorrelation lengths on P e (isotropic random field)

  Figure 6. Effect of L ln x on P e when L ln y =0.25

Table 1 . Shear Strength and Elastic Properties of Soil, Footing, and Interface

 1

	Variable	Soil	Footing Interface
	c	20kPa	N/A	20kPa
	φ	30 o	N/A	30 o
	ψ s =2/3 φ	20 o	N/A	20 o
	E	60MPa 25GPa	N/A
	ν	0.3	0.4	N/A
	K n	N/A	N/A	1GPa
	K s	N/A	N/A	1GPa

Table 2 . Evolution of the Failure Threshold with the Different Levels of the iSS and with the Number of Realizations N s (β=0, L ln x =2.5 and L ln y =0.25)

 2

	Failure threshold		Number of realizations at each level (N s )
	C j for each level j	200	400	600	800	1000	1200
	C 1	0.00191	0.00199 0.00189	0.00204	0.00191	0.00195
	C 2	0.00103	0.00099 0.00096	0.00110	0.00102	0.00103
	C 3	0.00041	0.00032 0.00021	0.00037	0.00036	0.00034
	C 4

Table 3 . Effect of the Radius of the Hypersphere on the Number of Realizations Required to Calculate P e (L ln x =2.5 and L ln y =0.25)

 3

		MCS	β=0 (Classical SS)	β=10	iSS β=11	β=11.5
	P e (x10 -4 )	3.40	3.65	3.58	3.36	3.45
	Number of levels	-	4	3	3	2
	number of realizations	30,000	3,700	2,800	2,800	1,900
	computation time (minutes) 210,000	25,900	19,600 19,600 13,300