I want to use Importance sampling to estimate probability of default of an insurance company within the next $t$ years. The company starts with capital $C$ at $t=0$. Each year it gains $p > 0$ in premiums and losses $X \sim F$ independently of other years:

$$Z_t = C + tp - \sum_{i=1}^{t} X_i$$

probability of default can be written as:

$$P(min_{\lt t \leq T}Z_t < 0) = E(1(min_{\lt t \leq T}Z_t < 0) )$$

Implemeting this in R , by letting x be a vector of simulated yearly losses. The function returns True if the company defaulted

default <- function(x,C,p){
  T <- length(x)
  Z_t<- C + p*seq(1,T,by=1) - cumsum(x)
  return(min(Z_t) < 0)


Let the loss distribution be independent $LogNormal(0,1/4)$. Simulating from this distribution and then taking the average However it is ineffeicient since the default probability very small ( less than $0.1\%$). When $t=20$ years, $C=15$ and $p=1$.

The idea is to instead sample from $Y_i \sim LogNormal(\mu,1/4)$. For some suitable choice of $\mu$.To make the default probability higher.

Using Importance sampling: $E(h(X)) = \int h(x)f_X(x)dx =\int h(x)f(x)\frac{f_Y(x)}{f_Y(x)}dx = E\Big(\frac{h(Y)f_X(Y)}{f_Y(Y)}\Big)$

we can take the expectation of: $$1(min_{\lt t \leq T}Z_t < 0) \frac{f_X(Y)}{f_Y(Y)} $$

The ratio of two lognormal densities $\frac{f_X(Y)}{f_Y(Y)}$ with $X \sim LogNormal(0,1/4), Y \sim LogNormal(\mu,1/4)$ with $\mu$ and the data as a parameters can be implemented as

w <- function(Y,mu){

Suppose now that we want to estimate Probability of default in the upcomming $t=20$ years with Starting capital $C=15$ and with premium $p=1$.

We run 10000 simulations with $\mu=0.5$ and finally average

N <- 10000
mu <- 0.5
           {y <- exp(rnorm(20,mean=mu,sd=1/2))

[1] 0.6254667

62% which is way too large. I cant see where this went wrong. Can someone give me some advice?


1 Answer 1


In the code

       y <- exp(rnorm(20,mean=mu,sd=1/2))

'y' is a vector of dimension 20. However,

w <- function(Y,mu){

returns a vector as well since the entry Y is vector of dimension 20. It should be a product

w <- function(Y,mu){

or equivalently

w <- function(Y,mu){

With this modification, the outcome is quite small

  > mean(hw1)
  [1] 0.000450381

as expected. And also highly variable as repeated calls to the code demonstrate.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.