Draft Conditioning by Permutation Monte Carlo for Continuous-Time - PowerPoint PPT Presentation

1 Draft Conditioning by Permutation Monte Carlo for Continuous-Time Markov Chains Pierre L’Ecuyer Universit´ e de Montr´ eal, Canada joint work with Zdravko Botev , New South Wales University, Australia Rohan Shah , The University of Queensland, Brisbane, Australia Bruno Tuffin , Inria–Rennes, France Eighth International Workshop on Simulation, Vienna, September 2015

2 Draft Object of this talk Think of a static system in which components have a random state, and there is a small probability u that their combined state is not good enough for the system to work properly. Estimating u by standard Monte Carlo can be very inefficient. In this talk we discuss a class of methods that transform the (simple) static model into a (more complicated) dynamic system, a continuous-time Markov chain, whose simulation can provide a much more efficient estimator of u when combined with conditional Monte Carlo.

3 Draft CTMC model Some rare-event estimation problems can be recast as follows. We have a continuous-time Markov chain (CTMC) { X ( γ ) , γ ≥ 0 } over state space X , with X (0) = x 0 ∈ X . In state x ∈ X , jump rate is λ ( x ) and next state is x ′ with (transition) probability p ( x ′ | x ). (Works also for a density p .)

3 Draft CTMC model Some rare-event estimation problems can be recast as follows. We have a continuous-time Markov chain (CTMC) { X ( γ ) , γ ≥ 0 } over state space X , with X (0) = x 0 ∈ X . In state x ∈ X , jump rate is λ ( x ) and next state is x ′ with (transition) probability p ( x ′ | x ). (Works also for a density p .) Let T 0 = 0 and T 1 < T 2 < · · · the jump times. At T j , chain jumps from state X j − 1 = X ( T j − 1 ) to state X j = X ( T j ). The associated discrete-time chain (DTMC) is { X j , j ≥ 0 } . Let ∆ ⊂ X (target) and C = inf { j ≥ 0 : X j ∈ ∆ } (a stopping time). We want to estimate some function of T C = inf { T j : X j ∈ ∆ } ; e.g., u ( γ ) = P [ T C > γ ]. Conditional on the DTMC trajectory H = ( X 0 , X 1 , . . . , X C ), the A j = T j − T j − 1 are independent exponential with rates Λ j = λ ( X j − 1 ).

4 Draft Monte Carlo (MC) for function of T C To estimate u = u ( γ ) = P [ T C > γ ], generate n realizations I 1 , . . . , I n of I = I [ T C > γ ], and compute estimator n I n = 1 ¯ � I i . n i =1 When u is very small ( { T C > γ } is a rare event), direct MC fails because the rare event occurs too rarely.

4 Draft Monte Carlo (MC) for function of T C To estimate u = u ( γ ) = P [ T C > γ ], generate n realizations I 1 , . . . , I n of I = I [ T C > γ ], and compute estimator n I n = 1 ¯ � I i . n i =1 When u is very small ( { T C > γ } is a rare event), direct MC fails because the rare event occurs too rarely. Relative error √ 1 − u MSE [¯ � I n ] RE [¯ I n ] def here = = √ nu → ∞ when u → 0 . u For example, if u ≈ 10 − 10 , we need n ≈ 10 12 to have RE [¯ I n ] ≤ 10%. We would like bounded RE (or almost) when u → 0.

5 Draft Conditioning on the DTMC Alternative: replace I = I [ T C > γ ] by W = P [ T C > γ | H ]. Idea used for static network reliability estimation by Elperin, Gertsbach, Lomonosov [1974, 1991, 1992, etc.]. Conditional on H , T C = A 1 + · · · + A C is hypoexponential, with cdf and pdf: F ( γ | H ) = P [ T C ≤ γ | H ] = 1 − e t 1 exp( Q γ ) 1 and f ( γ | H ) = − e t 1 exp( Q γ ) Q1 , where e t 1 = (1 , 0 , . . . , 0), 1 = (1 , . . . , 1) t , and   − Λ 1 Λ 1 0 · · · 0 . .   0 − Λ 2 Λ 2 0 .     ... ... Q = .   0 0 0    .  .   . 0 0 − Λ C − 1 Λ C − 1   0 · · · 0 0 − Λ C

6 Draft Higham (2009) has a stable and accurate algorithm for the matrix exponential that takes O ( C 3 ) time. Slow when C is large.

6 Draft Higham (2009) has a stable and accurate algorithm for the matrix exponential that takes O ( C 3 ) time. Slow when C is large. In some applications, Λ j decreases strictly with j . Then, explicit formulas ( O ( c 2 ) time): C C � � e − Λ j γ p j ( H ) = 1 − (1 − e − Λ j γ ) p j ( H ) W = P [ T C > γ | H ] = j =1 j =1 where C Λ k � p j ( H ) = . Λ k − Λ j k =1 , k � = j But does not work if some Λ k − Λ j = 0 (or near).

7 Draft Permutation Monte Carlo estimator To estimate u ( γ ) = P [ T C > γ ], sample n independent realizations of H and compute W = W ( H ) = P [ T C > γ | H ] for each, say W 1 , . . . , W n . PMC estimator: n W n = 1 ¯ � W i . n i =1 Since P [0 < W < 1] = 1, one has Var [ W ] < Var [ I ] and Var [ ¯ W n ] < Var [¯ I n ] (strictly) .

8 Draft Bounded relative error (BRE) Suppose the distribution of T C depends on a parameter ǫ and u = u ( ǫ ) = P [ T C ( ǫ ) > γ ] → 0 when ǫ → 0. Then, I n ] = Var [¯ I n ] = 1 − u RE 2 [¯ → ∞ when ǫ → 0 . u 2 nu

8 Draft Bounded relative error (BRE) Suppose the distribution of T C depends on a parameter ǫ and u = u ( ǫ ) = P [ T C ( ǫ ) > γ ] → 0 when ǫ → 0. Then, I n ] = Var [¯ I n ] = 1 − u RE 2 [¯ → ∞ when ǫ → 0 . u 2 nu W n ] = E [ W 2 / u 2 − 1] < ∞ iff E [ W 2 / u 2 ] ≤ K . However, RE 2 [ ¯ If the number of paths H is finite, this happens iff when ǫ → 0, P [ H ] W 2 ( H ) / u 2 ≤ K for all H . Holds iff P [ H ] �→ 0 for all H . A sufficient condition: whenever p ( x ′ | x ) > 0, p ( x ′ | x ) �→ 0 when ǫ > 0. More specific (weaker) conditions can be formulated on certain examples. Our development adapts easily to the case where we want to estimate P [ T C ( ǫ ) ≤ γ ] → 0 instead of P [ T C ( ǫ ) > γ ].

9 Draft Simple example: Credit risk An lending firm (or bank) has m obligors. With probability u i , obligor i defaults ( X i = 1) and the firm loses ℓ i from it. Total loss: m � L = L ( X ) = X i ℓ i . i =1 We want to estimate u = P [ L > ℓ ∗ ].

9 Draft Simple example: Credit risk An lending firm (or bank) has m obligors. With probability u i , obligor i defaults ( X i = 1) and the firm loses ℓ i from it. Total loss: m � L = L ( X ) = X i ℓ i . i =1 We want to estimate u = P [ L > ℓ ∗ ]. Exact formula (2 m terms): m � I [ L ( x ) > ℓ ∗ ] � u x i i (1 − u i ) 1 − x i . u = x ∈{ 0 , 1 } m i =1

9 Draft Simple example: Credit risk An lending firm (or bank) has m obligors. With probability u i , obligor i defaults ( X i = 1) and the firm loses ℓ i from it. Total loss: m � L = L ( X ) = X i ℓ i . i =1 We want to estimate u = P [ L > ℓ ∗ ]. Exact formula (2 m terms): m � I [ L ( x ) > ℓ ∗ ] � u x i i (1 − u i ) 1 − x i . u = x ∈{ 0 , 1 } m i =1 MC: generate X = ( X 1 , . . . , X m ), compute L and I = I [ L > ℓ ∗ ]. Repeat n times and take the average.

10 Draft From static model to CTMC PMC: define Y i exponential with rate λ i = − ln(1 − u i ). Then P [ Y i ≤ 1] = 1 − e − λ i = u i . Define a CTMC whose state is a binary vector X ( γ ) = ( X 1 ( γ ) , . . . , X m ( γ )) where X i (0) = 0 and X i ( γ ) jumps from 0 to 1 at time γ = Y i . We have P [ X i (1) = 1] = u i and P [ L ( X (1)) > ℓ ∗ ] = u .

10 Draft From static model to CTMC PMC: define Y i exponential with rate λ i = − ln(1 − u i ). Then P [ Y i ≤ 1] = 1 − e − λ i = u i . Define a CTMC whose state is a binary vector X ( γ ) = ( X 1 ( γ ) , . . . , X m ( γ )) where X i (0) = 0 and X i ( γ ) jumps from 0 to 1 at time γ = Y i . We have P [ X i (1) = 1] = u i and P [ L ( X (1)) > ℓ ∗ ] = u . Let π the permutation s.t. Y π (1) < · · · < Y π ( m ) (ordered jumps). Conditional on π , we can forget the Y i ’s, switch the X i ’s one by one until L ( X ) > ℓ ∗ , say at step C = min { j : L ( X j ) > ℓ ∗ } . The CTMC has initial jump rate Λ 1 = λ 1 + · · · + λ m . At the j th jump of the DTMC, X π ( j ) switches from 0 to 1 and the jump rate decreases to Λ j +1 = Λ j − λ π ( j ) . Thus, Λ j is strictly decreasing in j and C cannot exceed m . Easy to apply PMC.

11 Draft A CTMC that removes the defaults Instead of adding the defaults one by one, we can construct a CTMC that starts with all obligors in default, and remove the defaults one by one until L ( X ) ≤ ℓ ∗ , say at step C ′ . The “anti-default” that switches X i from 1 to 0 occurs at time R i , taken as exponential with rate µ i = − ln(1 − e − λ i ), so P [ R i ≤ 1] = P [ Y i > 1]. We have a similar PMC estimator of u , with C ′ , Λ ′ j , ..., defined from the µ i : W ′ = P [ A ′ 1 + · · · + A ′ C ′ > 1 | π ′ ] . In this example, we expect C ≪ C ′ , i.e., we should need much less than half the obligors to default, so the first approach should be faster. But the second approach may win in terms of variance.

12 Draft BRE, functional estimation, etc. Suppose u i = u i ( ǫ ) → 0 when ǫ → 0. If u i ( ǫ ) / u j ( ǫ ) ≥ δ > 0 for all i , j , then the PMC estimator has bounded relative error (sufficient condition for BRE).