# 1 - Assurance maritime

## Question 1

**A l’aide d’histogrammes et de graphes de probabilités, montrer que les lois normale et exponentielle ne sont pas des modèles plausibles pour modéliser les montants des sinistres.**

# Chargement du tableau de données
x <- scan("sinistres.csv")

# Tri des données (ordre croissant)
x <- sort(x)

# Nombre de données
n <- length(x)

# Histogramme à classes de même largeur
rangex <- max(x)-min(x)
a0 <- min(x) - 0.025*rangex
ak <- max(x) + 0.025*rangex

k <- round(log(n)/log(2)+1)
bornes <- seq(a0,ak,(ak-a0)/k)

hist(x, prob=T, breaks=bornes)

# Graphe de probabilités pour la loi normale
plot(x, qnorm(seq(1:n)/(n+1)))

# Graphe de probabilités pour la loi exponentielle
plot(x, log(1-seq(1:n)/(n+1)))

Dans le cas du graphe de probabilité de la loi normale, nous n'observons pas une droite.

Dans le cas du graphe de probabilité de la loi exponentielle, nous n'observons pas tout à fait à une droite. Une valeur exrême échappe à l'alignement.

## Question 2

# Fonction de répartition empirique
plot(ecdf(x))

# Indicateurs statistiques
summary(x)
var(x)

# Calcul des quantiles dans le cas d'une loi normale
pnorm(1.361, 2.168, 3.258) # 1er quantile
pnorm(2.322, 2.168, 3.258) # 3ème quantile

$$\bar{X_n} = 2.168, \bar{X_n}^{2} = 4.700, S_n'= 3.258$$

1er quartile : $1.361$
3ème quartile : $2.322$

Si $X_i \sim \cal{N}(2.168, 3.258)$, on devrait avoir $P(X_i \leq 1.361) = 0.25$ et $P(X_i \leq 2.322) = 0.75$.

Or $P(X_i \leq 1.361) = 0.402$ et $P(X_i \leq 2.322) = 0.519$.

Donc la loi normale n'est pas un bon modèle de par la distribution des données.

De plus, si la loi exponentielle était un bon modèle, on aurait $S_n' \approx \bar{X_n}^{2}$

## Question 3

On a la densité de la loi ${\cal P}a(a, b)$ : $$f(x) = \frac{a \, b^a}{x^{1+a}} \, \mathbb{1}_{[b,+\infty[}(x)$$

Soit $X$ une variable aléatoire de loi ${\cal P}a(a, b)$.

Sa fonction de répartition est la suivante : $$F_X(x) = \int_{-\infty}^{x} \frac{a \, b^a}{t^{1+a}} \, \mathbb{1}_{[b,+\infty[}(t)dt$$

Ainsi : $$F_X(x) = ab^{a}\int_{b}^{x} t^{-(1+a)}dt$$

Finalement, la fonction de répartition est : $$F_X(x) = ab^{a}\bigg(\frac{1}{ab^{a}} - \frac{1}{ax^{a}}\bigg) = 1 - \bigg(\frac{b}{x}\bigg)^{a}$$

Son espérance est : $$E(X) = \int_{-\infty}^{+\infty} tf(t)dt = ab^{a}\int_{b}^{+\infty} t^{-a}dt$$

On trouve donc : $$E(X) = ab^{a}\bigg[\frac{t^{1-a}}{1-a}\bigg]_{b}^{+\infty}$$

Finalement, l'espérance est : $$E(X) = \frac{ab^{a}x^{1-a}-ab}{1-a} = \frac{ab}{a-1}$$

De plus : $$E\big(X^{2}\big) = \frac{ab^{2}}{a-2}$$

Ainsi $$V(X) = E\big(X^{2}\big) - E(X)^{2} = \frac{ab^{2}}{a-2} - \bigg(\frac{ab}{a-1}\bigg)^{2} = \frac{ab^{2}}{(a-2)(a-1)^{2}}$$

Il est nécessaire que $a > 2$ pour que cette loi admette une espérance et une variance finies.

## Question 4

On pose $$ h : x \rightarrow ln(1 - x)$$
$$ g : x \rightarrow ln(x)$$
Alors : $$ h(F(x)) = -ag(x) + aln(b) $$
On note : $\alpha$ le coefficient directeur du graphe de probabilité et $\beta$ l'ordonnée à l'origine
Donc $$ a_g = -\alpha$$ et $$b_g = e^{\frac{\beta}{a_g}}$$

## Question 5

En divisant $V(X)$ par $E(X)^2$ qui est non nul on obtient :

$$ \frac{V(X)}{E(X)^2} (a^2-2a)= 1 $$

$$ \Leftrightarrow (a^2-2a) - \frac{E(X)^2}{V(X)} = 0 $$

Après résolution de l'équation du second degré on obtient :

Etant donné que $a > 2$ :

$$ a = 1 + \sqrt{1+\frac{E(X)^2}{V(X)}}$$

Et donc :

$$\tilde{a}_n = 1 + \sqrt{1+\frac{\bar{X}_n^2}{S_n'^2}}$$

Pour l'estimateur de b, on a :
$$ b = \frac{(a-1)E(X)}{a}$$
Donc en injectant l'expression de a en fonction de $E(X)$ et $V(X)$ dans cette dernière équation on obtient :
$$\tilde{b}_n = \frac{{\bar{X}_n} \sqrt{1+\frac{\bar{X}_n^2}{S_n'^2}}}{1 + \sqrt{1+\frac{\bar{X}_n^2}{S_n'^2}}}$$

## Question 6
On veut que $x_1,...,x_n \geq b$, on aura dans ce cas :
$$ {\cal L}(x_1,...,x_n,a,b) = \prod_{i=1}^n \frac{ab^a}{x_i^{1+a}}$$
sinon, on aurait :
$$ {\cal L}(x_1,...,x_n,a,b) = 0$$
${\cal L}$ est une fonction croissante en $b$ (car $a > 2$), donc pour maximiser ${\cal L}$ on prend $b$ le plus grand possible, c'est à dire on prend $b = min_{i}(x_i) = x_1^{*}$ (car on veut $x_1,...,x_n \geq b$).
Donc
$$\tilde{b}_n = min_{i}(x_i) = x_1^{*}$$
On a :
$$ ln({\cal L}(x_1,...,x_n,a,b)) = n(ln(a) + aln(b)) - (1+a) \sum_{i=1}^{n} ln(x_i)$$
Donc :
$$ \frac{\partial}{\partial a} ln({\cal L}(x_1,...,x_n,a,b)) = n\bigg(\frac{1}{a} + ln(b)\bigg) - \sum_{i=1}^{n} ln(x_i)$$
Donc, en égalisant cette dérivée partielle à 0, on obtient :

$$ \tilde{a}_n = \frac{n}{\sum_{i=1}^nln(\frac{x_i}{b})} $$

## Question 7

**La loi ${\cal P}a(a, b)$ vous parait-elle un modèle plausible pour les montants de sinistres ? Si oui, calculer toutes les estimations de $a$ et $b$ pour ces données. Le calcul théorique du biais et de la variance des estimateurs étant trop complexe, on déterminera quels sont les meilleurs estimateurs dans la partie 2, à l’aide de simulations.**

# définition de la fonction de répartition de la loi P_a(a,b)

repartition <- function(x, a, b) {
  return(1-((b/x)^a))
}

x <- scan("sinistres.csv")

X = log(x)
Y = log(1-repartition(x, 5, 7))

# Affiche le graphe de proba 
plot(X, Y)

Le graphe de probabilité nous donne une droite. La La loi ${\cal P}a(a, b)$ paraît être un modèle plausible.

x <- scan("sinistres.csv")
x <- sort(x)
xn <- mean(x)
vn <- var(x)
n <- length(x)

xg <- log(x)
yg <- log(1 - seq(1:n)/(n+1))
modele <- lm( yg ~ xg)
coefs <- coef(modele)

a4 <- -coefs[1]
b4 <- exp(coefs[2]/a4)

a5 <- 1 + sqrt(1 + (xn^2)/(vn))
b5 <- (xn*sqrt(1 + (xn^2)/(vn)))/(1 + sqrt(1 + (xn^2)/(vn)))

b6 <- min(x)
a6 <- n/sum(log(x/b6))

Nous trouvons $a_g = -0.379$, $b_g = 279$, $\tilde{a}_n = 2.563$, $\tilde{b}_n = 1.322$, $\bar{a}_n = 2.190$ et $\bar{b}_n = 1.201$.

# 2 - Vérifications expérimentales à base de simulations

## Question 1

On a $F_X(x) = P(X \leq x) = 1-\big(\frac{b}{x}\big)^{a}$.

$$P(Y \leq y) = P\bigg(ln\bigg(\frac{X}{b}\bigg) \leq y\bigg) = P\bigg(\frac{X}{b} \leq e^{y}\bigg) = P\bigg(X \leq be^{y}\bigg) = F_X(be^{y}) = 1-e^{-ay} = F_Y(y)$$

On peut donc simuler un échantillon de taille $n$ de la loi ${\cal P}a(a, b)$ à partir de la loi exponentielle $Exp(a)$.

<!-- Simuler les n échantillons à partir d'une loi exponentielle : rexp(n, rate = a) -->

## Question 2

a)

$$\frac{1}{a} = E(Y) $$

un estimateur de $\frac{1}{a}$ est :

$$\frac{1}{\tilde{a}_n}= \frac{1}{n} \sum_{i=1}^n ln\bigg(\frac{X}{b}\bigg)$$

on remarque que :

$$ \frac{n}{\tilde{a}_n} \sim \gamma(n, a) $$

et donc :

$$ \frac{2na}{\tilde{a}_n} \sim \gamma\bigg(n, \frac{1}{2}\bigg) \sim \chi^2_{2n} $$

or

$$ P\bigg(\frac{x\tilde{a}_n}{2n} \leq a \leq \frac{y\tilde{a}_n}{2n}\bigg) = P\bigg(x \leq \frac{2na}{\tilde{a}_n} \leq y\bigg) = F_{\chi^2_{2n}}(y)-F_{\chi^2_{2n}}(x) $$

on prend $x$ et $y$ tels que $F_{\chi^2_{2n}}(y)=1-\frac{\alpha}{2}$ et $F_{\chi^2_{2n}}(x)=\frac{\alpha}{2}$

$$y = z_{2n,\frac{\alpha}{2}}$$

$$x = z_{2n,1-\frac{\alpha}{2}}$$

Finalement

$$ P\bigg(\frac{z_{2n,1-\frac{\alpha}{2}}\tilde{a}_n}{2n} \leq a \leq \frac{z_{2n,\frac{\alpha}{2}}\tilde{a}_n}{2n}\bigg) = 1-\alpha $$

L'intervalle de confiance au seuil $\alpha$ est donc :

$$\bigg[\frac{z_{2n,1-\frac{\alpha}{2}}\tilde{a}_n}{2n},\frac{z_{2n,\frac{\alpha}{2}}\tilde{a}_n}{2n}\bigg]$$

b)

## Question 3

# {r include = False} pour ne pas afficher le code, mais uniquement les sorties

m <- 1000 # nombre d'échantillons
n <- 1000 # taille des échantillons
a <- 5
b <- 7

# estimateurs a et b graphiques (question 1.4)
a4 <- numeric(length=m)
b4 <- numeric(length=m)

# estimateurs a et b par la méthode des moments (question 1.5)
a5 <- numeric(length=m)
b5 <- numeric(length=m)

# estimateurs a et b par la méthode du maximum de vraisemblance (question 1.6)
a6 <- numeric(length=m)
b6 <- numeric(length=m)

for(i in 1:m)
{
    Y <- rexp(n, a)
    X <- b*exp(Y)
    xn <- mean(X) # moyenne empirique
    vn <- var(X) # variance éstimée
    xg <- log(sort(X))
    yg <- log(1 - seq(1:n)/(n+1))
    modele <- lm( yg ~ xg)
    coefs <- coef(modele)
    a4[i] <- -coefs[1]
    b4[i] <- exp(coefs[2]/a4[i])
    
    a5[i] <- 1 + sqrt(1 + (xn^2)/(vn))
    b5[i] <- (xn*sqrt(1 + (xn^2)/(vn)))/(1 + sqrt(1 + (xn^2)/(vn)))
    
    b6[i] <- min(X)
    a6[i] <- n/sum(log(X/b6[i]))
}

#calcul des biais de chaque estimateur
biais_a4 = mean(a4) - a
biais_a5 = mean(a5) - a
biais_a6 = mean(a6) - a
biais_b4 = mean(b4) - b
biais_b5 = mean(b5) - b
biais_b6 = mean(b6) - b

# calcul des erreurs quadratiques de chaque estimateur
err_quadra_a4 = var(a4) + (biais_a4)^2
err_quadra_b4 = var(b4) + (biais_b4)^2
err_quadra_a5 = var(a5) + (biais_a5)^2
err_quadra_b5 = var(b5) + (biais_b5)^2
err_quadra_a6 = var(a6) + (biais_a6)^2
err_quadra_b6 = var(b6) + (biais_b6)^2

# Enfin, on prend l'estimateur ayant le plus petit biais et la plus petite erreur quadratique moyenne pour trouver le meilleur

Les estimateurs calculés par la méthode du maximum de vraisemblance sont les meilleurs estimateurs de $a$ et de $b$. Leur biais et leur erreur quadratique sont les plus faibles.

## Question 4

a <- 5
b <- 7
epsilon <- 0.001
m <- 1000
#
prob_sup_a5 <- numeric(length = 91)
prob_sup_a6 <- numeric(length = 91)
for(n in seq(1000, 10000, by = 100))
{
    j = 1
    for(i in 1:m)
    {
        Y <- rexp(n, a)
        X <- b*exp(Y)
        xn <- mean(X) # moyenne empirique
        vn <- var(X) # variance éstimée
        a5 <- 1 + sqrt(1 + (xn^2)/(vn))
        a6 <- n/sum(log(X/min(X)))
        if(abs(a5 - a) > epsilon)
        {
            prob_sup_a5[j] <- prob_sup_a5[j] + 1
        }
        if(abs(a6 - a) > epsilon)
        {
            prob_sup_a6[j] <- prob_sup_a6[j] + 1
        }
        
    
    }
    prob_sup_a5[j] <- prob_sup_a5[j]/m
    prob_sup_a6[j] <- prob_sup_a6[j]/m
    j <- j + 1
}
plot(seq(1000, 10000, by = 100), prob_sup_a5)
plot(seq(1000, 10000, by = 100), prob_sup_a6)

## Question 5

# 

for(i in 5:1000)
{
    for(j in 1:1000)
    {
        Y <- rexp(i, a)
        X <- b*exp(Y)

        xn <- mean(X)
        vn <- var(X)

        a5 <- 1 + sqrt(1 + (xn^2)/(vn))
        b5 <- (xn*sqrt(1 + (xn^2)/(vn)))/(1 + sqrt(1 + (xn^2)/(vn)))

        b6 <- min(X)
        a6 <- n/sum(log(X/b6))
    }
}