Single-Index Importance Sampling with Stratification

Abstract

In many stochastic problems, the output of interest depends on an input random vector mainly through a single random variable (or index) via an appropriate univariate transformation of the input. We exploit this feature by proposing an importance sampling method that makes rare events more likely by changing the distribution of the chosen index. Further variance reduction is guaranteed by combining this single-index importance sampling approach with stratified sampling. The dimension-reduction effect of single-index importance sampling also enhances the effectiveness of quasi-Monte Carlo methods. The proposed method applies to a wide range of financial or risk management problems. We demonstrate its efficiency for estimating large loss probabilities of a credit portfolio under a normal and t-copula model and show that our method outperforms the current standard for these problems.

Appendix

1.1 Proofs

Proof of Proposition 1

The mean and variance follow from

$$\mathbb {E}_g(\hat{\mu }^{\tiny{\text{SIS}}}_n) = E_g(\Psi (\varvec{X}) w(T)) = \mathbb {E}_g(m(T)w(T)) = \int _{\Omega _g} m(t) \frac{f_T(t)}{g_T(t)}g_T(t)\;\mathrm {d}t=\mu _{\tiny{\text{SIS}}}$$

and

$$n {\text {Var}}_g(\hat{\mu }^{\tiny{\text{SIS}}}_n) + \mu _{\tiny{\text{SIS}}}^2 = \mathbb {E}_g\left( \Psi ^2(\varvec{X})w(T)\right) =\int _{\Omega _g} m^{(2)}(t)\frac{f_T^2(t)}{g_T^2(t)}\;\mathrm {d}t.$$

Asymptotic normality follows from the central limit theorem. Next, we need to find $g_T$ among all g that give unbiased estimators so that the variance, or equivalently $\mathbb {E}_g(m^{(2)}(T)w(T))$, is minimal when $\Psi (\varvec{x})\ge 0$ or $\Psi (\varvec{x})\le 0$ for all $\varvec{x}\in \Omega$. Let $\Omega _{\tiny{\text{ub}}}=\{t\in \Omega _f:m(t)f_T(t)\ne 0\}$. By Jensen’s inequality,

$$\begin{aligned} \mathbb {E}_g\left( m^{(2)}(T)w^2(T)\right)&\ge \left( \mathbb {E}_g\left( \sqrt{m^{(2)}(T)}w(T)\right) \right) ^2 \\&= \left( \int _{\Omega _g} \sqrt{m^{(2)}(t)}w(t)\;\mathrm {d}t\right) ^2 = \left( \int _{\Omega _f} \sqrt{m^{(2)}(t)}w(t)\;\mathrm {d}t\right) ^2 \end{aligned}$$

The last inequality follows since $\hat{\mu }^{\tiny{\text{SIS}}}_n$ is assumed to be unbiased, i.e., $\Omega _{\tiny{\text{ub}}}\subseteq \Omega _g$ and the fact that $\sqrt{m^{(2)}(t)}f_T(t)=0$ for $t\not \in \Omega _{\tiny{\text{ub}}}$ (as $m(t)=0$ implies $m^{(2)}(t)=0$ by the assumption on $\Psi$). The right hand side of the inequality is a constant independent of the choice of $g_T$, namely the minimum variance among all SIS estimators. To achieve equality, or equivalently to minimize the variance, set $g_T\propto \sqrt{m^{(2)}}(t)f_T(t)$ for $t\in \Omega _{\tiny{\text{ub}}}$ and the claim follows.

Proof of Proposition 2

Let $\Omega _T^{(i)}=\{t\in (t_{\inf },t_{\sup }): \lambda _{i} \le t < \lambda _{i+1}\}$ where $\lambda _i=G_T^\leftarrow ( (i+1)/n)$ and note that $\mathbb {P}(T\in \Omega _T^{(i)})=1/n$ for $i=1,\dots ,n$. Then

$$\begin{aligned} \mathbb {E}(\hat{\mu }^{\tiny{\text{SSIS}}}_n)&=\frac{1}{n}\sum _{i=1}^n \mathbb {E}_g\left( \Psi (\varvec{X})w(T)\mid T \in \Omega _T^{(i)}\right) =\frac{1}{n}\sum _{i=1}^n \mathbb {E}_g\left( \mathbb {E}_g(\Psi (\varvec{X})w(T)\mid T) \mid T \in \Omega _T^{(i)}\right) \\&= \frac{1}{n}\sum _{i=1}^n \int _{\lambda _i}^{\lambda _{i+1}}m(t) \frac{f_T(t)}{g_T(t)}g_T(t)\;\mathrm {d}t = \mu _{\tiny{\text{SIS}}} \end{aligned}$$

The expression for the variance is a slight generalization of (Glasserman et al. (1999), Lemma 4.1) in that stratification is combined with IS, bit it can be proved similarly. Let $\eta _{n}(t)$ denote the index i so that $t\in \Omega _T^{(i)}$. Then

$$\begin{aligned} n{\text {Var}}(\hat{\mu }^{\tiny{\text{SSIS}}}_n)=\frac{1}{n}\sum _{j=1}^n {\text {Var}}_g\left( \Psi (\varvec{X})w(T)\mid T \in \Omega _T^{(i)}\right) =\mathbb {E}_g\left( {\text {Var}}_g\left( \Psi (\varvec{X})w(T)\mid \eta _n(T)\right) \right) . \end{aligned}$$

Let $\xi =\mathbb {E}_g(\Psi (\varvec{X})w(T)\mid T)=m(T)w(T)$ and define the sequence $\xi _n=\mathbb {E}_g(\xi \mid \eta _n(T))$. Note that the $\sigma -$algebra generated by $\eta _n(T)$ forms an increasing family as n increases through a constant multiple of power two. Observe that $\mathbb {E}_g(|xi|)<\infty$ and $\sup _n \xi _n< \mathbb {E}_g(\Psi ^(\varvec{X})w^2(T))=\mathbb {E}_g(m^{(2)}(T)w^2(T))<\infty$. Also, $\xi _n$ is a martingale if n increases through a constant multiple of powers of two as it is a Doob’s martingale; (see Karlin and Taylor (1975), p. 246). Then using the arguments as in (Glasserman et al. (1999), Lemma 4.1), it follows that ${\text {Var}}_g(\hat{\mu }^{\tiny{\text{SSIS}}}_n)=\sigma _{\tiny{\text{SIS}}}^2/n+o(1)$.

The expression for the optimal density and variance expressions follow as in the proof of Prop. 1 by applying Jensen’s inequality. It remains to show that the SSIS estimator is asymptotically normal, which we show by applying the Lyapunov Central Theorem; (see Kole et al. (2007), p. 134). Let $m_i = \mathbb {E}_g(\Psi (\varvec{X})w(T)\mid T\in \Omega _T^{(i)})$ and $v_i^2={\text {Var}}_g(\Psi (\varvec{X})w(T)\mid T\in \Omega _T^{(i)})$. It is easily seen that $(1/n)\sum _{i=1}^n m_i=\mu_{\tiny{\text{SIS}}}$ and $(1/n)\sum _{i=1}^n v_i^2=\sigma _{\tiny{\text{SIS}}}^2+o(1)$. For any $i=1,\dots ,n$, we have

$$\begin{aligned}&\mathbb {E}_g\left( |\Psi (\varvec{X}_i)w(T_i)-m_i|^{2+\delta }\right) \le 2^{2+\delta }\left( \mathbb {E}_g\left( |\Psi (\varvec{X}_i)w(T_i)|^{2+\delta }\right) +\mathbb {E}_g\left( |m_i| ^{2+\delta }\right) \right) \\&= 2^{2+\delta } \left( \mathbb {E}_g\left( |\Psi (\varvec{X})w(T)|^{2+\delta } \mid T \in \Omega _T^{(i)}\right) + \mathbb {E}_g\left( |\mathbb {E}_g(\Psi (\varvec{X})w(T)\mid T\in \Omega _T^{(i)})|^{2+\delta }\right) \right) \\&\le 2^{2+\delta } \left( \mathbb {E}_g\left( |\Psi (\varvec{X})w(T)|^{2+\delta } \mid T \in \Omega _T^{(i)}\right) + \mathbb {E}_g\left( \mathbb {E}_g(|\Psi (\varvec{X})w(T)|^{2+\delta }\mid T\in \Omega _T^{(i)}) \right) \right) \\&=2^{3+\delta }\mathbb {E}_g\left( |\Psi (\varvec{X})w(T)|^{2+\delta }\mid T\in \Omega _T^{(i)}\right) , \end{aligned}$$

where the first inequality follows from the $c_{\tau }$ inequality as in (Loeve (1963), p. 155). The Lyapunov condition is satisfied, since

$$\begin{aligned}&\frac{1}{(\sum _{i=1}^n\sigma _i^2)^{1+\delta /2}} \sum _{i=1}^n\mathbb {E}_g\left( |\Psi (\varvec{X}_i)w(T_i)-m_i|^{2+\delta }\right) \\&\le \frac{2^{3+\delta }}{(\sum _{i=1}^n\sigma _i^2)^{1+\delta /2}} \sum _{i=1}^n\mathbb {E}_g\left( |\Psi (\varvec{X}_i)w(T_i)|^{2+\delta }\mid T\in \Omega _T^{(i)}\right) \\&= \frac{2^{3+\delta } n}{(n\sigma _{\tiny{\text{SSIS}}}^2+o(n))^{1+\delta /2}} \mathbb {E}_g\left( |\Psi (\varvec{X})w(T)|^{2+\delta }\right) \rightarrow 0,\quad n\rightarrow \infty , \end{aligned}$$

by the assumption. The Lyapunov Central Limit Theorem together with Slutsky’s Theorem implies $(\hat{\mu }^{\tiny{\text{SSIS}}}_n-\mu _{\tiny{\text{SIS}}})/\sqrt{n}\underset{}{\overset{\tiny {\text {d}}}{\rightarrow }}{\text {N}}(0,\sigma _{\tiny{\text{SSIS}}}^2)$.

Proof of Proposition 3

Recall that $T_i$ satisfies $T_i = G_T^\leftarrow ( (i+U_i-1)/n)$ where $U_i\overset{\tiny {\text {ind.}}}{\sim }\text {U}(0,1)$ for $i=1,\dots ,n$, and are therefore ordered, i.e., $T_1<T_2<\dots <T_n$. For any $i=1,\dots ,n$,

$$T_{i+1}-T_i=(G_T^{-1})'(\xi _i)\left( \frac{1+U_{i+1}-U_i}{n}\right) =\frac{1}{g_T(G_T^{-1}(\xi _i))}\left( \frac{1+U_{i+1}-U_i}{n}\right) =\mathcal {O}(1/n),$$

for some $\xi _i\in (T_i, T_{i+1})$, which implies that for any continuously differentiable function h, $h(T_{i+1})=h(T_i)+\mathcal {O}(1/n)$. Then we have

$$\begin{aligned} r_i^2&=\left( m(T_{i+1})+\varepsilon _{T_{i+1}} - m(T_i) - \varepsilon _{T_i}\right) ^2\\&= \left( m(T_{i+1})-m(T_i)\right) ^2 + \left( \varepsilon _{T_{i+1}}-\varepsilon _{T_i}\right) ^2-2(m(T_{i+1}) - m(T_i))(\varepsilon _{T_{i+1}}-\varepsilon _{T_i})\\&= (\varepsilon _{T_{i+1}}-\varepsilon _{T_i})^2 - 2(m(T_{i+1}) - m(T_i))(\varepsilon _{T_{i+1}}-\varepsilon _{T_i})+\mathcal {O}(1/n^2), \end{aligned}$$

and so

$$\begin{aligned} \mathbb {E}_g\left( r_i^2 w^2(T_i)\right)&= \mathbb {E}_g\left( \mathbb {E}_g(r_i^2w(T_i)\mid T_i, T_{i+1})\right) =\mathbb {E}_g\left( w^2(T_i)(v^2(T_i)+v^2(T_{i+1}))\right) +\mathcal {O}(1/n^2)\\&= 2 \mathbb {E}_g\left( w^2(T_i)v^2(T_i)\right) +\mathcal {O}(1/n), \end{aligned}$$

which means that

$$\begin{aligned} \mathbb {E}_g(\hat{\sigma }_{\tiny{\text{SSIS}}}^2)&= \frac{1}{2(n-1)}\sum _{i=1}^n \mathbb {E}_g(r_i^2w^2(T_i))=\frac{1}{n}\sum _{i=1}^n\mathbb {E}_g\left( v^2(T)w^2(T)\mid T\in \Omega _T^{(i)}\right) +\mathcal {O}(1/n)\\&= \mathbb {E}_g(v^2(T)w^2(T)) + \mathcal {O}(1/n) = \sigma _{\tiny{\text{SSIS}}}^2+ \mathcal {O}(1/n)\rightarrow \sigma _{\tiny{\text{SSIS}}}^2, \end{aligned}$$

which shows consistency.

Proof of Proposition 4

We use that $(\varvec{X}\mid T = t) \sim {\text {N}}_d(\varvec{\beta }t,I_d-\varvec{\beta }\varvec{\beta }^{\top })$(see Harris and Helvig (1965), Theorem 1) to compute the moment generating function of $\varvec{X}$. For $\varvec{a}\in \mathbb {R}^d$,

$$\begin{aligned} {\mathbb {E}_g(}&\mathbb {E}_g(\exp (\varvec{a}^{\top }\varvec{X})) = \mathbb {E}_g \left( \mathbb {E}(\exp (\varvec{a}^{\top }\varvec{X}) \mid T)\right) = \mathbb {E}_g\left( \exp \left( \varvec{a}^{\top }\varvec{\beta }T + \frac{1}{2} \varvec{a}^{\top }(I_d-\varvec{\beta }\varvec{\beta }^{\top })\varvec{a}\right) \right) \\&=\mathbb {E}_g(\exp (\varvec{a}^{\top }\varvec{\beta }T ))\exp \left( \frac{1}{2} \varvec{a}^{\top }(I_d-\varvec{\beta }\varvec{\beta }^{\top })\varvec{a})\right) =\exp \left( c\varvec{a}^{\top }\varvec{\beta }+\frac{1}{2}(\varvec{a}^{\top }\varvec{\beta })^2\sigma ^2\right) \times \\&\times \exp \left( \frac{1}{2} \varvec{a}^{\top }(I_d-\varvec{\beta }\varvec{\beta }^{\top })\varvec{a})\right) = \exp \left( \varvec{a}^{\top }(c\varvec{\beta }) + \frac{1}{2}\varvec{a}^{\top }(I_d+(\sigma ^2-1)\varvec{\beta }\varvec{\beta }^{\top }\varvec{a}\right) . \end{aligned}$$

By uniqueness of the moment generating function, $\varvec{X}\sim {\text {N}}_d(c\varvec{\beta }, I_d+(\sigma ^2-1)\varvec{\beta }\varvec{\beta }^{\top })$.