Category Archives: Course

Optimal Portfolios #1

This afternoon, I will start a crash course on financial portfolio optimization, with application in R. This week, we start with simple things, with the theoretical setup, without and with a risk free asset. We will discuss then the problem of estimating parameters, in a robust way. Then we introduce the idea of consider a more general criteria to quantify risk than the variance (but it means more general distributions… this point will be discussed further next time). The slides are available here, and R codes from there (in a Markdown)

Traffic Flow of Kota Kinabalu (with R)

This morning, we had our first practicals on network flows, using  an example mentioned in some papers published by Noraini Abdullah and Ting Kien Hua, max flow min cut theorem to minimize traffic congestion in Kota Kinabalu and application of the Shortest Path and Maximum Flow with Bottleneck in Traffic Flow of Kota Kinabalu. From the roads mentioned in the articles, I did try my best to locate the nodes on a map,

m=matrix(c(0,5.995910, 116.105520,
1,5.992737, 116.093718,
2,5.992066, 116.109883,
3,5.976947, 116.095760,
4,5.985766, 116.091580,
5,5.988940, 116.080112,
6,5.968318, 116.080764,
7,5.977454, 116.075460,
8,5.974226, 116.073604,
9,5.969651, 116.073753,
10,5.972341, 116.069270,
11,5.978818, 116.072880),3,12)

we can be visualized below

map = openmap(c(lat= 6.000, lon= 116.06),
c(lat= 5.960, lon= 116.12))
points(t(m[3:2,]),col="black", pch=19, cex=3 )

If the source is realistic (up north), I do not feel very confortable with the location of the sink (on the west). But let’s pretend it’s find (to do the maths, at least).

To extract information about edge capacity, on that network use the following code that will extract the three tables from the paper

location <- ''
out <- extract_tables(location)

with Windows, it seems to be necessary to download another package first

location <- ''
out <- extract_tables(location)

Now we can get out data frame with capacities[[2]])[[3]])

We can add those edges on our map (without the arrows to indicate the direction, it would be to heavy to read)

points(t(m[3:2,]),col="black", pch=19, cex=3 )
for(i in 1:nrow(E)){

To get the graph with capacities, an alternative is to use


but it does not respect geographical locations of nodes. It can actually be done using

plot(g, layout=as.matrix(B[,c("x","y")]))

To get a better understanding of the capacities of the road, use

plot(g, layout=as.matrix(B[,c("x","y")]),

From that network with capacities, the goal is to determine maximum flow on that network, from the source to the sink. This can be done with R using

> (m=max_flow(graph=g, source="s", target="t"))
[1] 2571

[1] 1191 1380 1422 1380 231 0 231 0 1149 1422 1149 0 0 1149 1422
[16] 1149

Our maximum flow is here 2571, which is different from was is actually claimed both in the two papers  max flow min cut theorem to… and application of the Shortest Path… (“the maximum flow for the capacitated network with 12 nodes and 16 edges of the selected scope in this study was 2598 vehicles per hour“) where there are clearly typos since values in the table and on the graph are different. Here I did use the ones from the tables.

plot(g, layout=as.matrix(B[,c("x","y")]),

That is nice, but rather odd. Actually, a much simpler flow can be considered, but the same global value

plot(g, layout=as.matrix(B[,c("x","y")]),

Nice, isn’t it. It is actually possible to do exactly the same on another paper they have, on the same city, traffic congestion problem of road networks in Kota Kinabalu.

location <- ''
out <- extract_tables(location)

Here the value of the maximal flow is 4017, just as they found in the original paper

Multinomial Logit as an Iterated Logit Regression

For the second section of the course at ENSAE, yesterday, we’ve seen how to run a multinomial logistic regression model. It is simply an extension of the binomial logistic regression. But actually, it is also possible to consider iterative binomial regressions.

Consider here a response variable Y with a multinomial distribution (3 factors to have something more general than the binomial), taking values \{A,B,C\}, with respective probabilities \mathbf{p}=(p_A,p_B,p_C). Here is a code to generate some multinomial variables

for(i in 1:B){Y[i]=sample(A,size=1,prob=C[i,])}

and here is a code to generate a dataset with n rows,


Let us generate a training dataset and a validation one


With a multivariate logistic regression

For convenience, consider the most popular factor in our training dataset

modalite=names(sort(table(DF1$Y),decreasing = TRUE))

Consider a regression model on the simulated dataset (with several covariates), let us estimate it, and let us get predictions.

reg=multinom(as.factor(Y) ~ ., data = DF1)
mp1=predict (reg, DF1, "probs")
mp2=predict (reg, DF2, "probs")

An alternative can be the following.
consider a first regression model on the Bernoulli variable Y_A=\mathbf{1}(Y=A). Actually, we will consider the most important factor, but for convenience, assume that it is A.
On our dataset, estimate that model, and get predictions. In the case where Y\neq A, define another Bernoulli variable Y_B=\mathbf{1}(Y=B|Y\neq A). We can estimate that model and derive two probabilities, \mathbb{P}(Y=B|Y\neq A) and \mathbb{P}(Y=C|Y\neq A) (the sum of the two being equal to 1). Based on those two models, it is possible to compute the three probabilities we are looking for. \mathbb{P}[Y=A] is obtained from the first model, and we can derive the other two from \mathbb{P}[Y=B|Y\neq A]\cdot\mathbb{P}[Y\neq A] and \mathbb{P}[Y=C|Y\neq A]\cdot\mathbb{P}[Y\neq A].

p11=predict (reg1, newdata=DF1, type="response")
p12=predict (reg2, newdata=DF1, type="response")
p21=predict (reg1, newdata=DF2, type="response")
p22=predict (reg2, newdata=DF2, type="response")

Let us compare the predicted probabilites, on the same dataset (here the training dataset)

> mmp1[1:9,c("0","1","2")]
0 1 2
1 0.19728737 0.4991805 0.3035321
2 0.17244580 0.5648537 0.2627005
3 0.19291753 0.5971058 0.2099767
4 0.09087176 0.7787304 0.1303978
5 0.23400225 0.4083022 0.3576955
6 0.18063647 0.6637352 0.1556283
7 0.13188881 0.7402710 0.1278401
8 0.13776970 0.6524959 0.2097344
9 0.12325864 0.6790336 0.1977078
> mp1[1:9,c("0","1","2")]
0 1 2
1 0.19691036 0.5022692 0.3008205
2 0.17123189 0.5680647 0.2607034
3 0.19293066 0.5984402 0.2086291
4 0.08821851 0.7813318 0.1304497
5 0.23470739 0.4109990 0.3542936
6 0.18249687 0.6602168 0.1572863
7 0.13128711 0.7400898 0.1286231
8 0.13525341 0.6553618 0.2093848
9 0.12090016 0.6815915 0.1975084

The two are very close. So yes, it is possible to see the multinomial regression as some sequential binomial regressions.

Nodal Regions and Flows

For practicals on networks and flows, we will use the R package flows dedicated to flows on networks

myflows <- prepflows(mat = nav, i = "i", j = "j", fij = "fij")
diag(myflows) <- 0

Select flows that represent at least 20% of the sum of outgoing flows for each urban area.

flowSel1 <- firstflows(mat = myflows/rowSums(myflows)*100, method = "xfirst",k = 20)

Then select the dominant flows (incoming flows criterion)

flowSel2 <- domflows(mat = myflows, w = colSums(myflows), k = 1)
flowSel <- myflows * flowSel1 * flowSel2
inflows <- data.frame(id = colnames(myflows), w = colSums(myflows))

and finally plot dominant flows map

opar <- par(mar = c(0,0,2,0))
sp::plot(GE, col = "#cceae7", border = NA)
plotMapDomFlows(mat = flowSel, spdf = UA, spdfid = "ID", w = inflows, wid = "id",wvar = "w", wcex = 0.05, add = TRUE,legend.flows.pos = "topright",legend.flows.title = "Nb. of commuters")
title("Dominant Flows of Commuters")

The code to get the background map is based on the GE object, defined in that package.

To go further on dominant flows read  Nystuen & Dacey (1961)

We will discuss in the last course, next week  two extensions that were not mentioned in the course. The first one is about congestion models. The second one is a nice application of flow to discuss sports issues in NBA (or NHL).

Traffic Flow of Kota Kinabalu (Malaysia)

For the second practicals of our course on networks and flows, we will study traffic flow of Kota Kinabalu (Malaysia), following several papers published by Noraini Abdullah and Ting Kien Hua, such as max flow min cut theorem to minimize traffic congestion in Kota Kinabalu, traffic congestion problem of road networks in Kota Kinabalu and application of the Shortest Path and Maximum Flow with Bottleneck in Traffic Flow of Kota Kinabalu.


Variables Catégorielles et Modèle Logistique

Petit complément, suite aux coquilles qu’il y avait dans les slides, sur une propriété de la régression logistique quand on régresse sur des variables catégorielles. On était sur la base des avocats

> avocat <- read.table("",header=TRUE,sep=",") > avocat$CLMSEX <- factor(avocat$CLMSEX, labels=c("M","F")) > avocat$MARITAL <- factor(avocat$MARITAL, labels=c("M","C","V","D")) > avocat=avocat[!$CLMSEX),]
> attach(avocat)
> sum((ATTORNEY==2)&(CLMSEX=="F"),na.rm=TRUE)/sum(CLMSEX=="F",na.rm=TRUE)
[1] 0.5256065
> sum((ATTORNEY==2)&(CLMSEX=="M"),na.rm=TRUE)/sum(CLMSEX=="M",na.rm=TRUE)
[1] 0.4453925

Autrement dit, la proportion d’hommes qui se sont fait représentés par un avocat est de 44.539%, et la proportion de femmes 52.56%. On peut visualiser le tableau croisé ci-dessous

> tab=xtabs(~ATTORNEY+CLMSEX,data=avocat)
> require(vcd)
> mosaic(tab, shade=TRUE, legend=TRUE)

Quand on fait une régression linéaire (Gaussienne), on retrouve ces probabilités dans les valeurs de coefficients: pour les femmes, on a 44.539%, et pour les hommes, on a la différence avec les femmes (cette dernière modalité étant la modalité de référence ici)

> reglm = lm((ATTORNEY==2) ~ CLMSEX, data=avocat)
> summary(reglm)
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.44539 0.02060 21.620 < 2e-16 ***
CLMSEXF     0.08021 0.02756 2.911  0.00367 **

C’est ce que l’on a ci-dessous

> sum(coefficients(reglm))
[1] 0.5256065
> coefficients(reglm)[1]

On a un résultat similaire avec une régression logistique

> reglogit = glm((ATTORNEY==2) ~ CLMSEX, data=avocat,family=binomial)
> summary(reglogit)
Estimate Std. Error z value Pr(>|z|)
(Intercept) -0.21930 0.08312 -2.639 0.00833 **
CLMSEXF      0.32182 0.11097  2.900 0.00373 **

même si les coefficients n’ont pas la même interprétation. En transformant ces coefficients, on retrouve très exactement les mêmes valeurs

> exp(sum(reglogit$coefficients[1:2])) /(1+exp(sum(reglogit$coefficients[1:2])))
[1] 0.5256065
> exp(reglogit$coefficients[1]) /(1+exp(reglogit$coefficients[1]))

Désolé pour la typo.

Actuariat de l’Assurance Non-Vie #3

Mardi, seconde partie du cours d’actuariat, avec les modèles de classification le matin, mais l’après midi, on devrait commencer les modèles de comptage. Les slides sont en ligne.

Le matin, je dois intervenir 15 minutes à l’IHP vers 9 heures, dans un colloque sur Artificial Intelligence for Fintech and Insurtech: si le RER traîne un peu trop, entre Luxembourg et Lozère, j’aurais peut-être 5 minutes de retard….

Articles for the Probability and Statistics Project

Here are some articles for the project for the graduate crash course on probability and statistics,

For those willing to work on datasets, consider success per school in a national exam, over time (here is the file)

base=read.table("", header=TRUE,sep=";",dec=",")

Some datasets will be uploaded soon

Métro: centralité et robustesse

Demain matin, nous aurons un TP pour le cours sur les réseaux, les flux et les transports. En particulier, en nous inspirant des travaux de Sybil Derrible, nous allons commencer par étudier la centralité dans les différents systèmes de métro, mais aussi la robustesse. Les matrices d’adjacence d’une trentaine de métros dans le monde sont en ligne dans un fichier xls. Histoire de gagner un peu de temps, le code pour créer une matrice d’adjacence peut être le suivant


On est ensuite prêt à manipuler le réseau,

iflo=graph_from_adjacency_matrix(Adj,mode = "undirected")

On va utiliser les notions vues en cours, sur la centralité, mais surtout, on travaillera sur Quantifying the robustness of metro networks, inspiré de The complexity and robustness of metro networks. Plusieurs fonctions utiles sont déjà programmées dans R, comme l’assortativité.



Base de Données, Devoir d’Actuariat

Pour les deux premiers devoirs du cours d’actuariat de l’assurance non-vie, on utilisera les mêmes données, dont le descriptif complet est en ligne ici. La première base d’apprentissage (sa version corrigée) est aussi en ligne

> base=read.csv2("",sep=",")
> dim(base)
[1] 100000 32
> table(base$survenance)

85636 14364

On reviendra en cours mardi prochain sur les attendus, mais l’idée sera de modélisation la variable survenance. Pour se second projet, on passera à la modélisation de la prime pure.

Actuariat de l’Assurance Non-Vie #2

Mardi, on commencera le cours d’actuariat, avec un cours d’introduction le matin, mais l’après midi, on attaquera plus sérieusement, avec une mise en bouche, avec les modèles de classification (de type {survenance} vs {non survenance} d’accident). Là encore, les slides sont en ligne.

Je parlerais aussi du premier devoir, qui portera sur les modèles que nous (re)verrons ici.


Actuariat de l’Assurance Non-Vie #1

Reprise du cours d’assurance non-vie à l’ENSAE la semaine prochaine. Les slides sont en ligne dans le fichier ENSAE-2017.7z, pour la globalité des slides, sinon je mettrais les liens au fur et à mesure (en corrigeant les typos, voire en rajoutant des sections, si besoin), avec cette semaine les slides pour l’introduction du cours. Et pour illustrer, j’ai emprunté quelques images au très beau livre Karambolage d’Arnold Odermatt.