actulab-2014-solution-cooperators/cooperators.Rnw

% 
\documentclass[letter]{article}
\usepackage[margin=0.75in]{geometry}
\usepackage{Sweave}
\usepackage{tikz} \usepackage{multirow}
\usepackage[utf8]{inputenc}
\usepackage{amsmath}
\numberwithin{equation}{section}
\usepackage{amsfonts} \usepackage{amsthm} \usepackage{thmtools}
\usepackage{hyperref}
\usepackage{cleveref}
\usepackage[hypcap]{caption}
\usepackage[off]{auto-pst-pdf}

\title{Risque de concentration en Alberta}
\author{François Pelletier}

\begin{document}

\maketitle

<<librairie,echo=FALSE>>=
set.seed(908144032)
library("actuar")
library("multicore")
source("mcsapply.R")
nb.conc <- 50
@

\section{Risque associé aux tempêtes de grêle}

Une tempête de grêle est un évènement de météo extrême qui peut créer beaucoup
de dommages dans une zone concentrée en peu de temps. Plusieurs facteurs
influencent la sévérité des dommages causés par ces tempêtes, dont le principal
est le diamètre moyen des grêlons. Le facteur influençant le nombre d'assurés
qui encourent un dommage est principalement la taille de l'orage de grêle et la
direction des vents.

J'ai choisi de modéliser les orages de grêle selon une force de 1 à 4, qui
influencera à la fois la taille moyenne de l'orage et la  taille des grêlons.
Afin de pouvoir obtenir facilement une formule analytique pour déterminer si un
assuré est touché par l'orage, j'ai choisi de lui donner une forme elliptique.
Ainsi, l'orage de force 4 sera formée de 4 ellipses concentriques, chacune ayant
un niveau de dommages allant de 1 à 4.

Je n'ai pas modélisé le lieu de formation des orages. Ainsi, un orage a la même
possibilité de se former à n'importe quel point à l'intérieur d'un rectangle
regroupant les assurés.

Les coordonnées en longitude et en latitude des assurés seront les seules
données utilisées ici.
La position des agents n'est pas nécessaire.
L'ensemble des montants et des paramètres des distributions statistiques ont été
attribués arbitraitement afin de donner des montants raisonnables à mon avis.

Enfin, afin de représenter le risque d'un assuré moins concentré, j'ai
simplement décidé de former un sous-groupe d'assurés formé d'un assuré sur
\Sexpr{nb.conc} dans l'ordre fourni, ce qui a pour effet de diminuer la concentration 
du risque tout en conservant l'effet de la densité de population, qui est déjà
concentrée sur le territoire.

\subsection{Définition des ellipses}

L'ellipse est une section conique qui peut être représentée selon une formule
paramétrique \ref{eq:ellipse} qui sert à définir les
coordonnées en fonction des paramêtres. L'ellipse est alors définie par son
centre $X_c,Y_c$ et la mesure de son plus grand rayon $a$ et de son plus petit
rayon $b$, qui sont perpendiculaires. L'inclinaison du plus grand rayon par
rapport à l'abscisse est donnée par l'angle $\varphi$ en radian.

\begin{eqnarray}
  \label{eq:ellipse}
  X(t)&=X_c + a\,\cos t\,\cos \varphi - b\,\sin t\,\sin\varphi\\
  Y(t)&=Y_c + a\,\cos t\,\sin \varphi + b\,\sin t\,\cos\varphi
\end{eqnarray}

La formule canonique \ref{eq:canonique} est une seconde forme qui sert à
déterminer si un point $(x,y)$ fait partie de l'ellipse.

\begin{eqnarray}
  \label{eq:canonique}
  \frac{\left((x-X_c)\cos \varphi-(y-Y_c)\sin \varphi\right)^2}{a^2} + 
  \frac{\left((x-X_c)\sin \varphi+(y-Y_c)\cos \varphi\right)^2}{b^2} = 1
\end{eqnarray}

Ces deux dernières équations nous permettrons de simuler des ellipses,
représentant les orages, et de déterminer si des assurés sont touchés par les
orages. Les functions suivantes seront utilisées à cette fin:

<<ellipse1>>=
source("ellipse.r")
@

\subsection{Importation des données}

On importe des données dans une structure
<<>>=
agents <- read.csv("agents.csv",sep=":")
assures <- read.csv("assures.csv",sep=":")
n.assures <- nrow(assures)
@

On visualise les données sur une carte. Les agents sont représentés par un
triangle rouge et les assurés par une croix noire.

<<fig=TRUE>>=
plot(assures$Longitude_Assure,assures$Latitude_Assure,pch=3,
		main="Carte des assurés et des agents",xlab="Longitude",
		ylab="Latitude")
points(agents$Longitude_Agent,agents$Latitude_Agent,pch=25,
		col="red")
@

On définit la zone où se situent les assurés par un rectangle:
<<>>=
coord_assures <- cbind(assures$Longitude_Assure,
		assures$Latitude_Assure)
range_longitude <- c(min(coord_assures[,1]),
		max(coord_assures[,1]))
range_latitude <- c(min(coord_assures[,2]),
		max(coord_assures[,2]))
@

On définit le groupe restreint pount l'évaluation avec moins de concentration
<<>>=
coord_assures2 <- coord_assures[seq(1,32669,by=nb.conc),]
@

La longitude est située entre \Sexpr{range_longitude[1]} et \Sexpr{range_longitude[2]}.
La latitude est située entre \Sexpr{range_latitude[1]} et \Sexpr{range_latitude[2]}.

\subsection{Hypothèses}

On pose plusieurs hypothèses de modélisation:

<<>>=
franchise <- 500
n.annees <- 100
n.moyen.orages <- 10
longueur.orage <- 1
largeur.orage <- .2
angle.dominant <- 1/2
ecart.angle <- 1
prob_forces <- c(.8, .15, .04, .01)
facteur_forces <- c(1,1.5,2,2.5)
param_dommages <- c(50,1000,20000,400000)
param_forme <- 2
@

\begin{enumerate}
  \item La franchise est de \Sexpr{franchise}\$.
  \item Le nombre d'années simulées est de \Sexpr{n.annees} années.
\item Le nombre d'orages par année suit une loi Poisson de moyenne
\Sexpr{n.moyen.orages}.
\item La longitude et la latitude du centre de l'orage sont distribués
uniformément sur le rectangle défini précédemment.
\item La longueur et la largeur suivent une loi exponentielle de paramètres
\Sexpr{1/longueur.orage} et \Sexpr{1/largeur.orage}. L'unité est en degrés de
latitude.
\item Les vents dominants sont de direction \Sexpr{angle.dominant} radians.
\item L'angle varie dans un invervalle de \Sexpr{ecart.angle} radian de part et
d'autre selon la distribution $\beta\left(2,2\right)$.
\item La force entre 1 et 4 suit une distribution discrète où les probabilités
sont respectivement \Sexpr{prob_forces[1]},
\Sexpr{prob_forces[2]}, \Sexpr{prob_forces[3]} et \Sexpr{prob_forces[4]}.
\item La force multiplie les dimensions de l'orage respectivement par
\Sexpr{facteur_forces[1]}, \Sexpr{facteur_forces[2]}, \Sexpr{facteur_forces[3]}
et \Sexpr{facteur_forces[4]}.
\item Les dommages suivent une distribution Pareto de paramètre de forme
\Sexpr{param_forme} et de paramètre d'échelle dépendant de la force, prenant les
valeurs
\Sexpr{param_dommages[1]}, \Sexpr{param_dommages[2]}, \Sexpr{param_dommages[3]}
et \Sexpr{param_dommages[4]}.
\end{enumerate}

\subsection{Simulation des orages}

En utilisant les hypothèse définies précédemment, on simule des orages

<<>>=
orages_annees <- rpois(n.annees,n.moyen.orages)
n.orages <- sum(orages_annees)
sim_xc <- runif(n.orages,range_longitude[1],range_longitude[2])
sim_yc <- runif(n.orages,range_latitude[1],range_latitude[2])
sim_a <- rexp(n.orages,1/longueur.orage)
sim_b <- rexp(n.orages,1/largeur.orage)
sim_angle <- angle.dominant - ecart.angle + 2*ecart.angle*
		rbeta(n.orages,2,2)
sim_force <- sample(x=1:4,size=n.orages,replace=TRUE,
		prob=prob_forces)
sim_param_orages_1 <- cbind(sim_xc,sim_yc,
		facteur_forces[1]*sim_a,
		facteur_forces[1]*sim_b,sim_angle)
sim_param_orages_2 <- cbind(sim_xc,sim_yc,
		facteur_forces[2]*sim_a,
		facteur_forces[2]*sim_b,sim_angle)
sim_param_orages_3 <- cbind(sim_xc,sim_yc,
		facteur_forces[3]*sim_a,
		facteur_forces[3]*sim_b,sim_angle)
sim_param_orages_4 <- cbind(sim_xc,sim_yc,
		facteur_forces[4]*sim_a,
		facteur_forces[4]*sim_b,sim_angle)
@

Pour chaque orage, on vérifie le nombre d'assurés touchés. On utilise une boucle
ici pour éviter de surcharger la mémoire vive. C'est plus rapide ainsi.

<<>>=
touches1 <- numeric(n.orages)
touches2 <- numeric(n.orages)
touches3 <- numeric(n.orages)
touches4 <- numeric(n.orages)
for(i in 1:n.orages)
{
	touches1[i] <- sum((sim_force[i] >= 1) * 
					dans_ellipse_param(coord_assures,
					sim_param_orages_1[i,]))
	touches2[i] <- sum((sim_force[i] >= 2) * 
					dans_ellipse_param(coord_assures,
					sim_param_orages_2[i,]))
	touches3[i] <- sum((sim_force[i] >= 3) * 
					dans_ellipse_param(coord_assures,
					sim_param_orages_3[i,]))
	touches4[i] <- sum((sim_force[i] >= 4) * 
					dans_ellipse_param(coord_assures,
					sim_param_orages_4[i,]))
}

rpareto_tronque <- function(n,shape,scale,deductible)
{
	pmax(rpareto(n,shape,scale)-deductible,0)
}

dommages <- mcsapply(mcsapply(as.list(touches1),
						rpareto_tronque,param_forme,
						param_dommages[1],franchise),sum)+
		mcsapply(mcsapply(as.list(touches2),
						rpareto_tronque,param_forme,
						param_dommages[2],franchise),sum)+
		mcsapply(mcsapply(as.list(touches3),
						rpareto_tronque,param_forme,
						param_dommages[3],franchise),sum)+
		mcsapply(mcsapply(as.list(touches4),
						rpareto_tronque,param_forme,
						param_dommages[4],franchise),sum)
dommages_cum <- cumsum(dommages)
orages_annees_cum <- cumsum(orages_annees)
dommages_annuels <- diff(dommages_cum[orages_annees_cum])
@

On obtient une moyenne de \Sexpr{mean(dommages_annuels)}, un écart-type de
\Sexpr{sd(dommages_annuels)} et les quantiles suivants:

<<>>=
quantile(dommages_annuels,c(.5,.75,.9,.95,.99))
@

\subsection{Comparaison avec moins de concentration}

On sélectionne un assuré sur \Sexpr{nb.conc} afin de réduire l'exposition et on
peut comparar les résultats avec les mêmes orages.

<<>>=
touches21 <- numeric(n.orages)
touches22 <- numeric(n.orages)
touches23 <- numeric(n.orages)
touches24 <- numeric(n.orages)
for(i in 1:n.orages)
{
	touches21[i] <- sum((sim_force[i] >= 1) * 
					dans_ellipse_param(coord_assures2,
					sim_param_orages_1[i,]))
	touches22[i] <- sum((sim_force[i] >= 2) * 
					dans_ellipse_param(coord_assures2,
					sim_param_orages_2[i,]))
	touches23[i] <- sum((sim_force[i] >= 3) * 
					dans_ellipse_param(coord_assures2,
					sim_param_orages_3[i,]))
	touches24[i] <- sum((sim_force[i] >= 4) * 
					dans_ellipse_param(coord_assures2,
					sim_param_orages_4[i,]))
}

dommages2 <- mcsapply(mcsapply(as.list(touches21),
						rpareto_tronque,param_forme,
						param_dommages[1],franchise),sum)+
		mcsapply(mcsapply(as.list(touches22),
						rpareto_tronque,param_forme,
						param_dommages[2],franchise),sum)+
		mcsapply(mcsapply(as.list(touches23),
						rpareto_tronque,param_forme,
						param_dommages[3],franchise),sum)+
		mcsapply(mcsapply(as.list(touches24),
						rpareto_tronque,param_forme,
						param_dommages[4],franchise),sum)
dommages2_cum <- cumsum(dommages2)
dommages2_annuels <- diff(dommages2_cum[orages_annees_cum])
@

On obtient une moyenne de \Sexpr{mean(dommages2_annuels)}, un écart-type de
\Sexpr{sd(dommages2_annuels)} et les quantiles suivants:

<<>>=
quantile(dommages2_annuels,c(.5,.75,.9,.95,.99))
@

\section{Risque associé au vol}

<<>>=
	taux_vol_petits <- 1600/100000
	taux_vol_grands <- 60/100000
@

\subsection{Hypothèses}
Le risque associé au vol peut être modélisé à l'aide d'une approche
fréquence-sévérité régulière. On pose donc comme hypothèse une distrubution
poisson pour la fréquence et une distribution gamma pour la sévérité.
On divise les vols en petits vols et en grands vols.
<<>>=
vols_petits <- taux_vol_petits * n.assures
freq_vols_petits <- rpois(n.annees,vols_petits)
vols_petits_totaux <- sum(freq_vols_petits)
sev_vols_petits <- rgamma(vols_petits_totaux,2,.005)
sev_vols_petits_cum <- cumsum(sev_vols_petits)
freq_petits_vols_cum <- cumsum(freq_vols_petits)
sev_vols_petits_annuels <- diff(sev_vols_petits_cum[freq_petits_vols_cum])

vols_grands <- taux_vol_grands * n.assures
freq_vols_grands <- rpois(n.annees,vols_grands)
vols_grands_totaux <- sum(freq_vols_grands)
sev_vols_grands <- rgamma(vols_grands_totaux,2,.0001)
sev_vols_grands_cum <- cumsum(sev_vols_grands)
freq_grands_vols_cum <- cumsum(freq_vols_grands)
sev_vols_grands_annuels <- diff(sev_vols_grands_cum[freq_grands_vols_cum])

vols_totaux <- rowSums(cbind(sev_vols_petits_annuels,sev_vols_grands_annuels))
@

On obtient une moyenne de vols annuels de \Sexpr{mean(vols_totaux)} avec un
écart-type de \Sexpr{sd(vols_totaux)}. Les quantiles de la distirubtion sont les
suivants:

<<>>=
quantile(vols_totaux,c(.5,.75,.9,.95,.99))
@



\end{document}