Freakonometrics

To content | To menu | To search

Miscellaneous › stochastic processes

Entries feed - Comments feed

Tuesday, March 27 2012

Ruin probability and infinite time

A couple of weeks ago, I had a discussion with a practitioner, working in some financial company, about ruin, and infinite time. And it remind me a weird result. Well, not a weird result, but a result I found disturbing, at first, when I was a student (that I rediscovered with the eyes of someone dealing with computational issues, seeing here a difficult theoretical question). Consider a simple ruin problem. A player has wealth . Then he flips a coin: tails he has a gain of 1, heads he experiences a loss of 1. At time , his wealth is where is associated to the th coin: is equal to 1 with probability (tails), and -1 with probability (heads). It is also possible to write

where can be interpreted as the net gain of the player. In order to get a good understanding of results that can be obtained. Assume to be given. Let denote the number of heads and the number of tails. Then , while . Let denote the number of paths to go from point A (wealth at time ) to point B (wealth at time ). Note that this is a Markovian problem, that can be modeled using Markov chains

But here, we will focus on combinatorial results. Hence,

In order to derive probabilities to reach , let denote the number of paths going from to . And let denote the number of paths going from to that do reach at some point between and . Using a simple reflexion property, then if and are positive,

Based on those reflexions, two results can be derived (focusing on probability, instead of counting paths). First, we can obtain that

(given that n and x have the same parity). The second result we can obtain is that

Based on those two expressions, if denotes the first time become null, given ,

then

This can be computed easily,

> x=10
> p=.55
> ProbN=function(n){
+ pb=0
+ if(abs(n-x) %% 2 == 0)
+ pb=x/n*choose(n,(n+x)/2)*(1-p)^((n+x)/2)*(p)^((n-x)/2)
+ return(pb)}
> plot(Vectorize(ProbN)(1:1000),type="s")

That looks nice... But if we look closer, we can wonder what

would be ? Since we have the distribution of a probabilty measure, we might expect one. But here

> sum(Vectorize(ProbN)(1:1000))
[1] 0.134385

And this is not due to calculation mistakes that we do not get 1 here. Actually, we should write

which might be interpreted as the probability of ruin, starting from , that we denote from now on. The term on the left can be approximated using monte-carlo simulations

> p=.55
> x=10
> m=1000
> simul=10000
> S=sample(c(-1,1),size=m*simul,replace=TRUE,prob=c(1-p,p))
> MS=matrix(S,simul,m)
> for(k in 2:m) MS[,k]=MS[,k]+MS[,k-1]
> T0=function(vm) which(vm<=(-x))[1]
> MTmin=apply(MS,1,T0)
> mean(is.na(MTmin)==FALSE)
[1] 0.1328

To check the validity of the relationship above, a simple (theoretical) recursive formula can be derived for the term on the right (ruin probability), namely

with a boundary conditions , and . Then is comes that

Note that it might be tricky to check using monte carlo simulation... since we cannot have an infinite number of runs. And we're dealing precisely with things that do occur when time is infinite. Actually, we can still check convergence, considering an upper limit for the number of runs, and then letting go to infinity. Note that an explicit formula can then be derived (using additional border condition )

Using the following code, it is possible to calculate ruin probability, in order to estimate .

> MSmin=apply(MS,1,min)
> mean(MSmin<=(-x))
[1] 0.1328
> (((1-p)/p)^x-((1-p)/p)^m)/(1-((1-p)/p)^m)
[1] 0.1344306

The following graph shows the evolution of ruin probability as a function of initial wealth (with monte carlo simulation, with a fixed horizon - including a confidence interval - versus the analytical expression)

Hence, with stopping times, one should remember that

and that those two terms can be approximated simply using simulations or standard approximations.

Wednesday, March 21 2012

La maudite constante des modèles ARIMA

Dans les modèles ARIMA, autant le dire tout de suite, les constantes c'est pénible ! Pour comprendre un peu mieux ce qui se passe, considérons ici un modèle ARMA avec constante. Ou pour commencer, juste un AR(1)
Supposons que la série  soit stationnaire, de moyenne . Alors en prenant l’espérance de part et d'autre dans l'expression précédente, , i.e. , ou . Plus généralement, avec un ARMA, , en prenant la encore l’espérance des deux cotés, on en déduit .

Simulons un processus AR(1),

Pour simuler des processus ARIMA, on pourrait utiliser la commande

> X=arima.sim(list(order=c(1,0,0),ar=1/3),n=1000)+2
> mean(X)
[1] 1.931767
mais comme le but est de comprendre ce qui se passe, autant faire les choses calmement,
> X=rep(NA,1010)
> X[1]=0
> for(t in 2:1010){X[t]=4/3+X[t-1]/3+rnorm(1)}
> X=X[-(1:10)]
> mean(X)
[1] 2.03397
I.e. ici le processus est de moyenne qui vaut 2.


Regardons maintenant ce que donnerait l'estimation du processus AR(1),

> arima(X, order = c(1, 0, 0))
 
Call:
arima(x = X, order = c(1, 0, 0))
 
Coefficients:
ar1  intercept
0.3738     2.0334
s.e.  0.0294     0.0487
 
sigma^2 estimated as 0.9318:  log likelihood = -1383.68
De manière un peu surprenante, le coefficient appelé intercept n'est pas la constante dans le modèle AR(1), mais la moyenne du processus. Autrement dit, R n'ajuste pas un processus
comme le laisserait penser l'intuition, mais un processus
Ces deux formes sont bien entendu équivalentes. Mais les coefficients estimés ne sont pas tout à fait ce que l'on attendait...
Plaçons nous maintenant dans le cas d'un processus non-stationnaire. L'extension naturelle serait de considérer un processus ARIMA(1,1,0), ou bien un processus tel que la différence soit un processus AR(1). Un processus ARIMA(1,1,0) avec une constante s'écrirait
en utilisant l'opérateur retard. Ceci fait penser à un modèle avec une tendance linéaire. Posons  (afin de se débarrasser de cette tendance). Alors
i.e.
ou encore
On note que  ou encore  suit alors un processus ARIMA(1,1,0) sans constante cette fois.
Afin de visualiser ce que donnerait l'inférence, simulons le processus suivant, qui est un processus AR(1) avec constante que l'on intègre:  avec
Peut-on retrouver les différents paramètres du processus avec R ?
Commençons (là encore) par simuler un tel processus,
> U=rep(NA,1010)
> U[1]=0
> for(t in 2:1010){U[t]=4/3+U[t-1]/3+rnorm(1)}
> U=U[-(1:10)]
> X=cumsum(U)
Sous R, on obtient l'estimation suivant si l'on tente de calibrer un modèle ARIMA(1,1,0)
> arima(X, order = c(1, 1, 0))
 
Call:
arima(x = X, order = c(1, 1, 0))
 
Coefficients:
ar1
0.8616
s.e.  0.0160
  
sigma^2 estimated as1.343:log likelihood = -1565.63

Ça ne convient pas du tout.... On peut tenter un processus AR(1) (avec constante) sur la série différenciée...

> arima(diff(X), order = c(1, 0, 0))
 
Call:
arima(x = diff(X), order = c(1, 0, 0))
 
Coefficients:
ar1  intercept
0.3564     2.0200
s.e.  0.0295     0.0486
 
sigma^2 estimated as 0.9782:  log likelihood = -1406.6
On progresse, sauf que comme auparavant, le terme qui est donne n'est pas la constante dans le modèle ARIMA, mais la moyenne du processus différencié. Mais cette fois, on a un interprétation, c'est que la constante est la pente de la tendance ! Si on estime la pente associée a , on récupère la même valeur:
> arima(X, order = c(1, 1, 0), xreg=1:length(X))
 
Call:
arima(x = X, order = c(1, 1, 0), xreg = 1:length(X))
 
Coefficients:
ar1  1:length(X)
0.3566       2.0519
s.e.  0.0296       0.0487
 
sigma^2 estimated as 0.9787:  log likelihood = -1406.82
Sur la figure ci-dessous, on retrouve le fait qu'en enlevant la tendance linéaire à  donneune série intégrée, sans constante (qui fait penser à une marche aléatoire)

Autrement dit

  • dans le cas d'une série stationnaire, la constante estimée n'est pas du tout la constante, mais la moyenne de la série temporelle
  • dans le cas d'une série non-stationnaire, la constante estimée dans la série différenciée a du sens, au sens ou il s'agit de la pente de la tendance (linéaire) du processus
Avant de conclure, une petite remarque. Quid de la prévision ? Si on commencer par une révision à l'aide du premier processus ARIMA, on obtient une prédiction pour avec un énorme intervalle de confiance (sans aucun bon sens)
> ARIMA1=arima(X, order = c(1, 1, 0))
> ARIMA2=arima(X, order = c(1, 1, 0), xreg=1:length(X))
> Xp1=predict(ARIMA1,20)
> Xp2=predict(ARIMA2,20,newxreg=
+ (length(X)+1):(length(X)+20))
> plot(960:1000,X[960:1000],xlim=c(960,1020),type="l")
> polygon(c(1001:1020,rev(1001:1020)),
+ c(Xp1$pred+2*Xp1$se,rev(Xp1$pred-2*Xp1$se)),
+ col=CL[3],border=NA)
> lines(1001:1020,Xp1$pred,col="red",lwd=2)
intervalle qui se visualise sur le graphique ci-dessous,

Si on regarde pour l'autre modèle,
> lines(1001:1020,Xp2$pred,col="blue",lwd=2)

Tous ceux qui auront reconnu ici des problèmes qui se posent dans la modélisation de la série http://freakonometrics.blog.free.fr/public/maths/viekt.png dans le modèle de Lee & Carter (1992) auront compris qu'il existe une vraie application a tout ce que je viens de raconter (je pense au particulier au commentaire poste par @ClaudeT en début de semaine) car la série des http://freakonometrics.blog.free.fr/public/maths/viekt.png ressemble très fortement à ce genre de série, non-stationnaire avec une tendance linéaire. Et qu'un modèle mal spécifié laisse à penser que l'incertitude est beaucoup beaucoup plus grande que ce qu'elle est vraiment (en plus de donner une prédiction surprenante à long terme !). Donc avant de faire tourner rapidement les codes, il vaut mieux regarder attentivement ce qu'ils font réellement...

Monday, March 19 2012

Simulation d'un processus de Lévy, et discrétisation

Avec @renaudjf, on discutait l'autre jour de la simulation d'un processus de Lévy. Et on se posait la question d'un algorithme optimal pour combiner un processus de Poisson (ou un process Poisson composé) avec un processus de Wiener (avec éventuellement un drift, voire une diffusion plus générale). En fait, pour générer des processus de Poisson, j'ai toujours eu l'habitude de simuler les durées entre sauts (avec des lois exponentielles, indépendantes, comme dans des vieux billets). Jean François me suggérait d'utiliser une propriété d'uniformité des sauts sur un intervalle de temps donné, conditionnellement aux nombres de sauts.

Commençons par la première piste. On peut générer un processus de Wiener, éventuellement avec un drift, et à coté, on peut générer les lois exponentielles  (qui vont correspondre aux durées entre sauts), et éventuellement aussi des amplitudes de sauts (e.g. des pertes qui suivent des lois exponentielles). On a ici

où . On commence par générer  en notant que

où les incréments  sont Gaussiens (centrés et de variance ) et indépendants les uns des autres. Quant aux durées entre sauts, les  , ce sont des lois exponentielles indépendantes, de moyenne . Voilà le code qui permet de générer les trois composantes,

n=1000
h=1/n
lambda=5
set.seed(2)
W=c(0,cumsum(rnorm(n,sd=sqrt(h))))
W=rexp(100,lambda)
N=sum(cumsum(W)<1)
T=cumsum(W[1:N])
X=-rexp(N)

Le hic est que pour le processus de Wiener, on a du discrétiser, alors que pour le processus de Poisson composé, non. Pourtant, il va bien se ramener sur une même échelle de temps. Une première piste est de créer vraiment la fonction 

Lt=function(t){
W[trunc(n*t)+1]+sum(X[T<=t])+lambda*t
}

et pour faire un dessin ensuite, c'est un jeu d'enfant

L=Vectorize(Lt)
u=seq(0,1,length=n+1)
plot(u,L(u),type="l",col="blue")

Une autre possibilité est d'utiliser une propriété d'uniformité du processus de Poisson que j'évoquais en introduction. Car le processus de Poisson vérifie une propriété remarquable: si  est la date où survient le ième saut, , alors conditionnellement au fait que , les variables  correspondent aux statistiques d'ordres de  variables indépendantes, uniformément distribuées sur , i.e.

Cette propriété se trouve dans Wolff (1982). L'idée de la démonstration est relativement simple. Commençons avec un (unique) saut, alors pour ,

i.e. on retrouve la fonction de répartition d'une loi uniforme sur . On itère ensuite avec 2 sauts, 3 sauts, etc.

La traduction en R de cette idée est tout simplement (car on se place sur )

N=rpois(1,lambda)
T=runif(N)
X=-rexp(N)

Ensuite, une stratégie est de discrétiser le processus de Poisson, avec le même pas de temps que le processus de Wiener,

indice=trunc(T*n)+1
saut=rep(0,n+1)
saut[indice]=X
processus=W+cumsum(saut)+lambda*u

On retrouve la même trajectoire qu'auparavant

plot(u,processus,type="l",col="red")

,,

Sauf qu'on a eu de la chance. Avec cette procédure, il ne faut pas que l'on ait deux sauts dans le même intervalle de temps ! Bon, il est vrai qu'une caractérisation du processus de Poisson est que

et donc on doit avoir très peu de chance d'avoir deux sauts au même instant d'autant plus que le pas de temps est petit. Mais "peu de chance" ne veut pas dire nul, et si on génère des milliers de trajectoires, la probabilité d'avoir une fois un soucis n'est pas négligeable.

Jean-François a eu l'idée brillante de proposer de tirer non pas des lois uniformes sur , mais des lois uniformes discrètes, dans

 .

T=c(0,sort(sample((1:(n-1)/n),size=N,replace=FALSE)))

sans remise afin d'éviter d'avoir deux sauts au même moment. L'idée était séduisante, mais il me restait à etre convaincu que les durées entre sauts... continuaient à être (presque) exponentielles.

Pour ça, on peut faire quelques tests. Par exemple, générer quelques simulations pour avoir une centaines de sauts (et donc une centaine de durées entre sauts), puis faire un test de loi exponentielle (de moyenne )

VT=0
for(ns in 1:20){
N=rpois(1,lambda)
if(N>0){
T=c(0,sort(sample((1:(n-1)/n),size=N,replace=FALSE)))
VT=c(VT,diff(T))
}}

On fait ici 20 boucles car on avait fixé

lambda=5

et j'avais dit que je voulais une centaine d'observations pour faire un test de loi (ce qui est purement arbitraire, on en conviendra). On peut ensuite faire un test d'ajustement de loi exponentielle, 

ks.test(VT[-1],"pexp",lambda)$p.value

Si on répète un grand nombre de fois, en changeant le pas de temps (ou le nombre de subdivisons de l'intervalle de temps), on notera qu'effectivement, avec un grand pas de temps (à gauche ci-dessous) on rejette souvent, voire presque toujours, l'hypothèse de loi exponentielle. Mais qu'assez vite, c'est une hypothèse qui n'est pas invraisemblable,

Je ne sais pas entre ces deux fonctions ce qui serait le plus rapide, mais on a deux jolis algorithmes pour générer les processus de Lévy. Non ?