# Estimate default by importance sampling (using R)

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){
exp(-log(Y)^2/0.5)/exp(-((log(Y)-mu)^2)/0.5)
}


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
hw1<-replicate(N,
{y <- exp(rnorm(20,mean=mu,sd=1/2))
default(y,C=15,p=1)*w(y,mu)}
)
mean(hw1)

[1] 0.6254667


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

In the code

       y <- exp(rnorm(20,mean=mu,sd=1/2))
default(y,C=15,p=1)*w(y,mu)


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

w <- function(Y,mu){
exp(-log(Y)^2/0.5)/exp(-((log(Y)-mu)^2)/0.5)}


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

w <- function(Y,mu){
prod(exp(-(log(Y)^2-(log(Y)-mu)^2)/0.5))}


or equivalently

w <- function(Y,mu){
exp(sum((log(Y)-mu)^2-log(Y)^2))/0.5)}


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.