© Revista de Matemática: Teoría y Aplicaciones 2022 29(2) : 239–260 cimpa – ucr issn: 1409-2433 (Print), 2215-3373 (Online)

doi: https:// doi. org/ 10.15517/ rmta. v29i2. 47938

estimación máxima verosimilitud de la probabilidad de ruina en el modelo de

riesgo clásico con reclamaciones exponenciales

maximum likelihood estimation of ruin probability in the classical risk model

with exponential claims

Ernesto A. Guerrero-Lara^{∗}  Jesús E. López-Flores^{†}

Henry G. Pantí-Trejo^{‡}

Received: 10/Sep/2021; Revised: 31/May/2022; Accepted: 22/Jun/2022

∗Universidad Autónoma de Yucatán, Facultad de Matemáticas, Mérida, México. E-Mail: ernesto.guerrero@correo.uady.mx

†Universidad Autónoma de Yucatán, Facultad de Matemáticas, Mérida, México. E-Mail: enrique_lopz.f@hotmail.com

‡Universidad Autónoma de Yucatán, Facultad de Matemáticas, Mérida, México. E-Mail: henry.panti@correo.uady.mx

240 e.a. guerrero-lara — j.e. lópez-flores — h.g. pantí-trejo

Resumen

Se calculan los estimadores de máxima verosimilitud para los parámetros que definen al proceso de Poisson compuesto en el proceso de riesgo clásico con reclamaciones exponenciales. Se prueba consistencia y normalidad asintótica de los estimadores obtenidos. Finalmente, con ayuda de la propiedad de invarianza de los estimadores de máxima verosimilitud, la normalidad asintótica y el método delta, se realiza una estimación puntual y por intervalos de la probabilidad de ruina.

Palabras clave: estimación máxima verosimilitud; probabilidad de ruina; modelo clásico de ruina; método delta.

Abstract

Maximum likelihood estimators are calculated for the parameters that define the compound Poisson process in the classical risk process with exponential claims. It is proved consistency and asymptotic normality for estimators obtained. Finally, with the help of invariance property of the maximum likelihood estimators, asymptotic normality and delta method, point and interval estimation of the ruin probability is performed.

Keywords: ruin probability; maximum likelihood estimation; classical ruin model; delta method.

Mathematics Subject Classification: 62F12, 91B05, 62M05, 62P05.

1 Introducción

estimación máxima verosimilitud de la probabilidad de ... 241

Las aseguradoras son empresas o instituciones que se encargan de proteger a un particular o un bien contra los riesgos que se puedan presentar a cambio del pago de una prima. Las primas representan la mayor parte de los ingresos de las aseguradoras y en la mayoría de los ramos de no vida, su cálculo está en función de la Siniestralidad Agregada de la cartera. Debido al comportamiento aleatorio presente en la siniestralidad agregada, es posible que en algún momento ésta supere al monto de capital acumulado con el que cuenta la aseguradora dando origen a lo que se conoce desde el punto de vista teórico como ruina. Por tal razón, es de interés el evento ocurre ruina y el cálculo de su probabilidad. Para ello se consideran modelos matemáticos que modelan el capital de la aseguradora a través del tiempo, los cuales son conocidos como modelos estocásticos. El modelo estocástico ampliamente utilizado para modelar la siniestralidad agregada es el proceso de Poisson compuesto que considera que el número de reclamaciones en cierto tiempo se distribuye como una variable aleatoria Poisson y que los montos de las reclamaciones son variables aleatorias independientes e idénticamente distribuidas entre sí e independientes del número de reclamaciones.

El primer trabajo conocido que propone un modelo de Poisson compuesto para modelar la siniestralidad fue por parte de Filip Lundberg en el año 1903 [7], básicamente consideraba a la siniestralidad agregada en cualquier intervalo de tiempo como un proceso de Poisson compuesto. Posteriormente en 1930 Harald Cramér, presentó una extensión a la propuesta de Lundberg, integrando su investigación sobre la aproximación de distribuciones de sumas de variables aleatorias independientes, trabajos posteriores contribuyeron de manera constante en esta área (ver [2]). En la actualidad diversas adaptaciones y generalizaciones del modelo antes mencionado han surgido (ver [1], [6], [12] y referencias que ahí aparecen).

El modelo basado en la teoría desarrollada por Lundberg se le conoce como proceso de riesgo clásico o modelo de Cramér-Lundberg. De acuerdo con [5], el proceso de riesgo clásico es una simplificación de la realidad, no obstante, es un modelo útil que puede proporcionar una perspectiva de las características de la operación del seguro. El caso más conocido y ampliamente estudiado de este modelo es cuando los montos de las reclamaciones siguen una distribución exponencial ([1], [6], [10], [12]). En este caso son dos los parámetros que definen al proceso de riesgo clásico: la intensidad del proceso de Poisson y la media de los montos de las reclamaciones.

El objetivo de este trabajo es presentar de una forma detallada el procedimiento de estimación puntual, vía máxima verosimilitud, de los parámetros del proceso de riesgo clásico, bajo el supuesto de reclamaciones exponenciales. De igual manera se demuestran propiedades de consistencia y normalidad asintótica de los estimadores obtenidos, siendo esta última, desde nuestro punto de vista, una demostración diferente a las que se pudieran encontrar en la literatura y una de las principales aportaciones de este trabajo. Además, el tener de manera explícita la probabilidad de ruina en el caso con reclamaciones exponenciales, permite ilustrar, de una forma más clara, el procedimiento de estimación puntual y, a través del método delta, una estimación asintótica por intervalos para la probabilidad de ruina. Finalmente, se proporciona un ejemplo numérico de los resultados obtenidos, se calcula el error cuadrático medio de las estimaciones y los intervalos de confianza asintóticos para la probabilidad de ruina. Lo anterior complementa lo realizado en [9] para el caso exponencial, el cual no fue estudiado a detalle en ese trabajo.

El artículo está organizado como sigue. En la Sección 2 se introduce el modelo clásico de Cramér-Lundberg y las cantidades que son relevantes para abordar este tema. En la Sección 3 se establecen los estimadores de máxima verosimilitud para el proceso que define al modelo clásico, se prueba la consistencia y normalidad asintótica de los estimadores obtenidos en la sección anterior. Posteriormente, se establece un estimador puntual para la probabilidad de ruina y se calculan, utilizando el método delta, intervalos asintóticos de confianza para la probabilidad de ruina. Finalmente en la Sección 4, se presenta un ejemplo numérico con los resultados previamente establecidos.

2 Modelo clásico de Cramér-Lundberg

Un proceso de Poisson N = {N_{t}, t ≥ 0} con intensidad λ > 0, es un proceso con incrementos independientes y estacionarios que satisface N_{0} = 0 y

P (N_{t} = n) = e^{−λt} \frac{(λt)^{n}}{n!} , n = 0, 1, …, t > 0.

Si además se considera una sucesión de variables aleatorias {Y_{n}} independientes e idénticamente distribuidas e independiente del proceso de Poisson N, se define el proceso de Poisson compuesto {S_{t}, t ≥ 0} con intensidad λ como

S_{t} = \begin{cases}0, & si N_{t} = 0; \\ ∑_{k=1}^{N_{t}} Y_{k}, & si N_{t} > 0. \end{cases}

El proceso de Poisson compuesto también posee incrementos independientes y estacionarios. Además, si E[Y_{i}] = μ, se satisface para toda t > 0

E(N_{t}) = λt, E (S_{t}) = λμt.  (1)

El proceso de riesgo clásico, también conocido como modelo clásico de Cramér-Lundberg, modela la evolución del capital de una aseguradora a lo largo del tiempo, considerando un capital inicial u, un ingreso de primas constante a una tasa c por unidad de tiempo y una siniestralidad agregada que se comporta de acuerdo a un proceso de Poisson compuesto con intensidad λ, donde la sucesión de variables aleatorias {Y_{n}} se asumen no negativas. De manera precisa, el modelo de Cramér-Lundberg U = {U(t, u), t ≥ 0} está dado por

U(t, u) = u + ct − S_{t}, t ≥ 0.  (2)

En la Figura 1 se presenta una posible trayectoria del modelo de Cramér-Lundberg, en el que los decrementos se presentan en el momento de que se presente una reclamación y el tamaño del decremento es precisamente el monto de la reclamación.

Figura 1: Trayectoria del modelo de Cramér-Lundberg.

Uno de los eventos de interés para este proceso es el tiempo en el cual la compañía aseguradora no cuenta con la solvencia para afrontar sus riesgos, esto es, cuando el capital de la compañía aseguradora es negativa, en términos del modelo de Cramér-Lundberg la ruina ocurre si el capital de la aseguradora U(t, u) es menor o igual a cero, por lo que el tiempo de ruina, el cual es un tiempo aleatorio, está dado por

τ(u) = min{t > 0 : U(t, u) ≤ 0}.

A su vez, la probabilidad de ruina, denotada por ψ(u), se define como

ψ(u) = P (τ(u) < ∞).

Un supuesto importante en el modelo de Cramér-Lundberg es la condición de beneficio neto, que indica que la cantidad de prima a cobrar c debe ser mayor a la prima de riesgo, es decir, el ingreso por primas debe ser mayor al valor esperado de los siniestros que se deben pagar, esto es, c > λμ. Si esta condición no se satisface entonces ψ(u) = 1, para toda u ≥ 0 [10]. En la Figura 1 esta condición se refleja por la tendencia creciente de la trayectoria.

Una de las cantidades que permiten aproximar la probabilidad de ruina o en su caso, calcularla de forma explícita, es el coeficiente de ajuste. De manera precisa, el coeficiente de ajuste es la solución positiva r^{∗} (si existe) de la ecuación

λM_{Y} (r) − λ − cr = 0,  (3)

donde M_{Y} denota la función generadora de momentos común de las variables de la sucesión {Y_{n}}, siempre que ésta exista. El coeficiente de ajuste r^{∗} se relaciona con la probabilidad de ruina de las siguientes dos maneras:

• Desigualdad de Lundberg: Para toda u ≥ 0,

ψ(u) ≤ e^{−r^{∗}u}.

• Aproximación de Crámer-Lundberg:

ψ(u) ∼ \frac{c − λμ}{λM′_{Y} (r^{∗}) − c}e^{−r^{∗}u}, cuando u → ∞,

donde M′_{Y} es la derivada de la función generadora de momentos común de las variables de la sucesión {Y_{n}} (ver [13], [5], [12]).

En general el coeficiente de ajuste no siempre existe, se tienen condiciones que permiten asegurar la existencia del coeficiente de ajuste r^{∗}. Por ejemplo en [5] se establece que si existe γ, con 0 < γ ≤ ∞, tal que M_{Y} (r) es finita para toda r < γ y lim_{r→γ} M_{Y} (r) = ∞, entonces existe una única solución positiva de la

ecuación (3). De manera particular, si {Y_{n}} tienen una distribución exponencial con media μ, el coeficiente de ajuste existe y resolviendo la ecuación (3) se verifica que

r^{∗} = \frac{1}{μ} − \frac{λ}{c}.

Por otro lado, no es una labor sencilla y en ocasiones no es posible obtener una fórmula explícita de la probabilidad de ruina. Existen diversas formas para calcular la probabilidad de ruina, desde expresiones que involucran series hasta ecuaciones integro-diferenciales. Para el caso de reclamaciones exponenciales con media μ, es posible proporcionar una expresión explícita para la probabilidad de ruina (ver [1], [10], [12]). De manera precisa, si los montos de las reclamaciones siguen una distribución exponencial con media μ entonces

ψ(u) = \frac{λμ}{c}e^{−(1/μ−λ/c)u}, u ≥ 0.  (4)

Para concluir esta sección, es importante notar que en el modelo de Cramér-Lundberg son dos los parámetros desconocidos que determinan su comportamiento probabilístico, a saber, λ y μ. Por lo que de ahora en adelante escribiremos θ = (λ, μ), para denotar al vector de parámetros a estimar.

3 Reclamaciones con distribución exponencial

Supongamos que la sucesión de los montos de las reclamaciones {Y_{n}, n ≥ 1}, tienen densidad exponencial con media μ:

f(y; μ) = \frac{1}{μ}e^{−y/μ}, y > 0.  (5)

Fijemos un tiempo T > 0 y observemos en el intervalo [0, T] al proceso de riesgo U definido en (2). Tenemos que N_{T} es el número total de reclamaciones en el intervalo [0, T] y Y_{1}, …, Y_{N_{T}} son los montos correspondientes de las reclamaciones individuales realizadas en el mismo intervalo. De esta manera, la muestra está dada por Y_{1}, …, Y_{N_{T}} , N_{T} y sus respectivos valores los denotaremos por y_{1}, …, y_{n}, n.

Es posible que no hayan reclamaciones en el intervalo [0, T], esto puede significar que la intensidad del proceso de Poisson λ es pequeña. Con el fin de establecer un procedimiento válido de estimación para los parámetros λ y μ, a lo largo de este artículo supondremos que se observa al menos una reclamación en el intervalo [0, T].

3.1 Estimación de parámetros

La estimación de θ = (λ, μ) se realiza por el método de máxima verosimilitud por lo que se procede a calcular la función de verosimilitud de la muestra. De la definición de proceso de Poisson compuesto se sigue que condicional a N_{T} = n, las variables aleatorias Y_{1}, …, Y_{N_{T}} tienen la misma distribución que Y_{1}, …, Y_{n}, y son independientes e idénticamente distribuidas. Considerando que la función de densidad de las Y_{i} está dada por la ecuacón (5), se tiene que la función de verosimilitud es:

L(λ, μ) = L(λ, μ; y_{1}, …, y_{n}, n)

= f_{Y_{1},...,Y_{N_{T}} ,N_{T}} (y_{1}, …, y_{n}, n)

= f_{Y_{1},...,Y_{N_{T}} |N_{T} =n}(y_{1}, …, y_{n}|N_{T} = n)P (N_{T} = n)

= (∏_{i=1}^{n} \frac{1}{μ}e^{−y_{i}/μ}) e^{−λT} \frac{(λT)^{n}}{n!}

= \frac{1}{μ^{n}} exp \{−\frac{1}{μ} ∑_{i=1}^{n} y_{i}\} e^{−λT} \frac{(λT)^{n}}{n!} .

Aplicando logaritmo se obtiene la función de logverosimilitud, la cual se puede escribir de la manera siguiente

ℓ(λ, μ) = ℓ_{1}(λ) + ℓ_{2}(μ) + cte,

donde

ℓ_{1}(λ) = −λT + n log λ
ℓ_{2}(μ) = −n log μ − \frac{1}{μ} ∑_{i=1}^{n} y_{i}
cte = − log n! + n log T

De aquí, las derivadas parciales de ℓ están dadas por

\frac{∂}{∂λ}ℓ(λ, μ) = \frac{∂}{∂λ}ℓ_{1}(λ) = −T + \frac{n}{λ}

\frac{∂}{∂μ}ℓ(λ, μ) = \frac{∂}{∂μ}ℓ_{2}(μ) = −\frac{n}{μ} + \frac{1}{μ^{2}} ∑_{i=1}^{n} y_{i}

De lo anterior notamos que \frac{∂^{2}}{∂λ∂μ}ℓ(λ, μ) = \frac{∂^{2}}{∂μ∂λ}ℓ(λ, μ) = 0 y las segundas

derivadas parciales satisfacen

\begin{matrix}\frac{∂^{2}}{∂λ^{2}} ℓ(λ, μ) = −\frac{n}{λ^{2}} \\ \frac{∂^{2}}{∂μ^{2}} ℓ(λ, μ) = \frac{n}{μ^{2}} − \frac{2}{μ^{3}} ∑_{i=1}^{n} y_{i} \end{matrix}  (6)

Resolviendo el sistema

\frac{∂}{∂λ}ℓ(λ, μ) = 0

\frac{∂}{∂μ}ℓ(λ, μ) = 0

y con ayuda de las ecuaciones (6), se puede verificar que λ̂, μ̂ dados por

λ̂ = \frac{N_{T}}{T} , μ̂ = \frac{1}{N_{T}} ∑_{i=1}^{N_{T}} Y_{i}  (7)

son los Estimadores de Máxima Verosimilitud (EMV) para λ, μ, respectivamente (ver por ejemplo [8] para criterios de maximización de funciones bivariadas).

Con ayuda de las ecuaciones en (1) se puede verificar que aunque los EMV para λ y μ consideran el comportamiento probabilístico de la muestra, coinciden con los que se obtienen por el método de momentos que sólo consideran medidas resumidas de ésta. Lo anterior proporciona un mayor sustento teórico a los estimadores obtenidos para λ y μ en esta sección.

3.2 Comportamiento asintótico de los EMV

En esta sección se establecen propiedades asintóticas para los EMV de λ y μ. Para tal fin, primero recordemos algunos conceptos sobre convergencia de sucesiones de vectores aleatorios.

Sea X_{1}, X_{2}, … una sucesión de vectores aleatorios de dimensión k. Decimos que la sucesión {X_{n}} converge a X

• casi seguramente (c.s.), denotado por X_{n} ⟶^{c.s.} X, si P (X_{n} → X) = 1;

• en probabilidad, denotado por X_{n} ⟶^{p} X, si para toda ε > 0, se satisface

lim_{n→∞} P (|X_{n} − X| < ε) = 1;

• en distribución, denotado por X_{n} ⟶^{d} X, si para cada punto de continuidad de la función de distribución conjunta F de X, F_{n}(x) → F (x), donde F_{n} son las funciones de distribución conjunta de X_{n}.

Estos tipos de convergencia están relacionados de la siguiente manera: la convergencia casi segura implica la convergencia en probabilidad y ésta a su vez implica la convergencia en distribución.

Ahora bien, si θ̂_{1}, θ̂_{2}, … es una sucesión de estimadores para el vector de parámetros θ satisfaciendo θ̂_{n} →^{c.s.} θ entonces se dice que la sucesión θ̂_{n} es consistente fuerte. Por otro lado, si se cumple θ̂_{n} →^{p} θ entonces la sucesión θ̂_{n} es consistente débil. En el caso \sqrt{n}(θ̂_{n} − θ) →^{d} N_{k}(0, Σ), donde N_{k}(0, Σ) denota un vector aleatorio normal multivariado de dimensión k, con vector de medias 0 y matriz de covarianzas Σ, se dice que la sucesión de estimadores θ̂_{n} es asintóticamente normal. La consistencia fuerte (débil) permite asegurar que, conforme el tamaño de muestra crece, las estimaciones se aproximan con probabilidad uno (con probabilidad alta) al verdadero valor del parámetro, mientras que la propiedad de asintoticidad normal permite aproximar la distribución (muestral) de los estimadores a través de una distribución normal.

Es conocido que cuando T → ∞,

\frac{N_{T}}{T} ⟶^{c.s} λ.

Por otro lado, la convergencia N_{T} ⟶^{c.s} ∞ y la Ley Fuerte de los Grandes Números implican

\frac{1}{N_{T}} ∑_{i=1}^{N_{T}} Y_{i} ⟶^{c.s} μ,

cuando T → ∞. Lo anterior permite establecer la siguiente proposición.

Proposición 3.1 Los EMV de λ y μ dados por:

λ̂_{T} = \frac{N_{T}}{T} , μ̂_{T} = \frac{1}{N_{T}} ∑_{i=1}^{N_{T}} Y_{i},

son consistentes.

Aunque la propiedad de normalidad asintótica de los EMV es bastante conocida, en aras de exhaustividad, se proporciona una demostración alterna de este hecho para los EMV de λ y μ, es decir, se proporciona una demostración de que

\sqrt{T} (λ̂_{T} − λ, μ̂_{T} − μ)^{t} ⟶^{d} N_{2}(0, I_{1}^{−1}(θ)),  (8)

donde I_{1}^{−1}(θ) es la inversa de la matriz de información de Fisher. Utilizando la definición de matriz de información de Fisher dada en [4]:

I(θ) = E [−\frac{∂^{2}}{∂θ_{i}∂θ_{j}} ℓ(θ)] ,

donde θ es el vector de parámetros del modelo, se calcula esta matriz para el modelo de Cramér-Lundberg con reclamaciones distribuidas exponencialmente

I_{T} (θ) = (\begin{matrix}E [−\frac{∂^{2}}{∂λ^{2}} ℓ(λ, μ)] & E [−\frac{∂^{2}}{∂λ∂μ}ℓ(λ, μ)] \\ E [−\frac{∂^{2}}{∂μ∂λ}ℓ(λ, μ)] & E [−\frac{∂^{2}}{∂μ^{2}} ℓ(λ, μ)] \end{matrix}) ,

donde θ = (λ, μ). Utilizando las derivadas calculadas en (6) y las identidades en (1), se verifica

I_{T} (θ) = T diag(λ^{−1}, λ/μ^{2}).

De aquí,

I_{1}^{−1}(θ) = diag(λ, μ^{2}/λ).  (9)

Así, la convergencia en (8) es equivalente a la siguiente convergencia en distribución

I1^{\frac{1}{2}} (θ)\sqrt{T} (λ̂_{T} − λ, μ̂_{T} − μ)^{t} ⟶^{d} N_{2}(0, I).  (10)

Sean

X_{T}^{1} = \frac{N_{T} − λT}{\sqrt{λT}} , X_{T}^{2} = \frac{\overline{Y}_{T} − μ}{μ/\sqrt{λT}},  (11)

donde

\overline{Y}_{T} = \frac{1}{N_{T}} ∑_{i=1}^{N_{T}} Y_{i}.

Por lo tanto, la convergencia en distribución (10), en términos del vector aleatorio (X_{T}^{1} , X_{T}^{2} )^{t} se escribe como:

(X_{T}^{1} , X_{T}^{2} )^{t} ⟶^{d} N_{2}(0, I).  (12)

Es conocido que cuando T → ∞, X_{T}^{1} ⟶^{d} Z y X_{T}^{2} ⟶^{d} Z, donde Z es una variable aleatoria normal estándar. Es decir se tienen las convergencias en distribución unidimensionales, sin embargo, esto no asegura la convergencia en distribución conjunta a un vector aleatorio normal bivariado estándar. El Teorema 3.3 asegura esto último y para demostrarlo, primero se demuestra el siguiente lema.

Lema 3.2 Sea X_{T}^{2} definido en (11). Entonces

E (e^{iθ_{2}X_{T}^{2}} | N_{T} ) ⟶^{c.s.} e^{− \frac{1}{2} θ_{2}^{2}},

cuando T → ∞.

Demostración. Sean

Y_{i}^{∗} = \frac{Y_{i} − μ}{μ}, i = 1, 2, ….

De esta manera,

X_{T}^{2} = \frac{\sqrt{λT}}{N_{T}} ∑_{i=1}^{N_{T}} Y_{i}^{∗}.

Ya que {Y_{n}^{∗}, n ≥ 1} también es una sucesión de variables aleatorias independientes, idénticamente distribuidas e independiente del proceso de Poisson N, se obtiene

E (e^{iθ_{2}X_{T}^{2}} | N_{T} = n) = (φ_{Y_{i}^{∗}} (\frac{θ_{2}\sqrt{λT}}{n} ))^{n} = (φ_{Y_{i}^{∗}} (\frac{θ_{2}α_{n}}{\sqrt{n}} ))^{n} ,

donde φ_{Y_{i}^{∗}} es la función característica de la variable Y_{i}^{∗} y α_{n} = \sqrt{\frac{λT}{n}}.

Ahora bien, las convergencias N_{T} →^{c.s.} ∞ y \frac{N_{T}}{T} →^{c.s.} λ, cuando T → ∞, aseguran que si T → ∞, entonces n → ∞, α_{n}^{2} → 1 y \frac{α_{n}^{2}}{n} → 0.

Por otro lado, ya que Y_{i} es una variable aleatoria exponencial, E(Y_{i}^{∗}) = 0,

E(Y_{i}^{∗2}) = 1. Por lo anterior y que \frac{α_{n}^{2}}{n} → 0, se sigue

φ_{Y_{i}^{∗}} (\frac{θ_{2}α_{n}}{\sqrt{n}} ) = 1 − \frac{θ_{2}^{2}α_{n}^{2}}{2n} + o (θ_{2}^{2} \frac{α_{n}^{2}}{n} ) .

De esta manera,

E (e^{iθ_{2}X_{T}^{2}} | N_{T} = n) = (1 − \frac{θ_{2}^{2}α_{n}^{2}}{2n} + o (θ_{2}^{2} \frac{α_{n}^{2}}{n} ))^{n} .

De la expresión anterior y debido a que N_{T} →^{c.s.} ∞, α_{N_{T}}^{2} →^{c.s.} 1, \frac{α_{N_{T}}^{2}}{N_{T}} →^{c.s.} 0, cuando T → ∞, se concluye

E (e^{iθ_{2}X_{T}^{2}} | N_{T} ) = (1 − (\frac{θ_{2}^{2}α_{N_{T}}^{2}}{2} ) \frac{1}{N_{T}} + o (θ_{2}^{2} \frac{α_{N_{T}}^{2}}{N_{T}} ))^{N_{T}} ⟶^{c.s.} e^{− \frac{1}{2} θ_{2}^{2}},

cuando T → ∞.

Ahora se está listo para demostrar el siguiente teorema.

Teorema 3.3 Sean

X_{T}^{1} = \frac{N_{T} − λT}{\sqrt{λT}} , X_{T}^{2} = \frac{\overline{Y}_{T} − μ}{μ/\sqrt{λT}},

donde

\overline{Y}_{T} = \frac{1}{N_{T}} ∑_{i=1}^{N_{T}} Y_{i}.

Entonces, cuando T → ∞,

(X_{T}^{1} , X_{T}^{2} )^{t} ⟶^{d} N_{2}(0, I).

Demostración. Se tiene que

E (e^{iθ_{1}X_{T}^{1} +iθ_{2}X_{T}^{2}} | N_{T} ) = e^{iθ_{1}X_{T}^{1}} E (e^{iθ_{2}X_{T}^{2}} | N_{T} ) .

Por otro lado, cuando T → ∞,

X_{T}^{1} ⟶^{d} Z,

donde Z es una variable aleatoria normal estándar. Este hecho, el Lema 3.2 y el Teorema de Slutsky aseguran cuando T → ∞ que

E (e^{iθ_{1}X_{T}^{1} +iθ_{2}X_{T}^{2}} | N_{T} ) = e^{iθ_{1}X_{T}^{1}} E (e^{iθ_{2}X_{T}^{2}} | N_{T} ) ⟶^{d} e^{iθ_{1}Z}e^{− \frac{1}{2} θ_{2}^{2}}.

Por lo tanto, si T → ∞

E (e^{iθ_{1}X_{T}^{1} +iθ_{2}X_{T}^{2}} ) = E (E (e^{iθ_{1}X_{T}^{1} +iθ_{2}X_{T}^{2}} | N_{T} ))
→ E (e^{iθ_{1}Z}e^{− \frac{1}{2} θ_{2}^{2}})
= e^{− \frac{1}{2} (θ_{1}^{2}+θ_{2}^{2})}.

Lo que muestra el resultado.

Del teorema anterior y de lo discutido previamente, se obtiene el siguiente corolario.

Corolario 3.4 Los EMV de λ y μ son asintóticamente normales. De manera precisa, se satisface

\sqrt{T} (λ̂_{T} − λ, μ̂_{T} − μ)^{t} ⟶^{d} N_{2}(0, I_{1}^{−1}(θ)),

donde

I_{1}^{−1}(θ) = diag(λ, μ^{2}/λ).

3.3 Estimación puntual y por intervalos de la probabilidad de ruina

En [3, Sección 7.2] se establece que si (θ̂_{1}, …, θ̂_{k}) es EMV de (θ_{1}, …, θ_{k}) y si τ(θ_{1}, …, θ_{k}) es cualquier función de los parámetros, entonces el EMV de τ(θ_{1}, …, θ_{k}) es τ(θ̂_{1}, …, θ̂_{k}). Esto permite asegurar que si λ̂ y μ̂ son los EMV de λ y μ, respectivamente, entonces el EMV de la probabilidad de ruina para reclamaciones exponenciales, tomando como función la ecuación (4), está dada por

ψ̂(u) = \frac{λ̂μ̂}{c}e^{−(1/μ̂−λ̂/c)u}.  (13)

Ya que la probabilidad de ruina es una función de los parámetros y su estimación una función de los EMV, los cuales son asintóticamente normales, es posible establecer, gracias al Método Delta, una fórmula para el cálculo de un intervalo asintótico de confianza para la probabilidad de ruina. Con fines ilustrativos, se enuncia el Método Delta Multivariado (ver [3], [11], [14]):

Teorema 3.5 (Método Delta Multivariado) Sea X = (X_{1}, …, X_{k}) un vector de variables aleatorias que satisface X ⟶^{p} Θ, donde Θ = (θ_{1}, …, θ_{k}). Además

\sqrt{n}(X − Θ) ⟶^{d} N_{k}(0, Σ),

siendo Σ la matriz de covarianzas del vector X. Entonces para una función g : ℝ^{k} → ℝ, con primeras derivadas parciales continuas con respecto a Θ se cumple

\sqrt{n}(g(X) − g(Θ)) ⟶^{d} N_{1}(0, g′(Θ)^{t}Σg′(Θ)),

donde g′(Θ) = (\frac{∂}{∂θ_{1}} g(Θ), …, \frac{∂}{∂θ_{k}} g(Θ))^{t}.

Por la Proposición 3.1 los EMV para θ = (λ, μ) son consistentes, por el Corolario 3.4 son normales asintóticamente, por la ecuación (4) para un capital inicial fijo u la probabilidad de ruina ψ(u) es una función de θ = (λ, μ) y además

\frac{∂}{∂λ}ψ(u) = \frac{μ}{c} (1 + \frac{λu}{c}) e^{−(1/μ−λ/c)u}, \frac{∂}{∂μ}ψ(u) = \frac{λ}{c} (1 + \frac{u}{μ}) e^{−(1/μ−λ/c)u}.

De esta forma se satisfacen las condiciones del Teorema 3.5, lo que implica

\sqrt{T} (ψ̂(u) − ψ(u)) ⟶^{d} N_{1}(0, ψ′(u)^{t}I_{1}^{−1}(θ)ψ′(u)), T → ∞,  (14)

donde ψ′(u) = (\frac{∂}{∂λ}ψ(u), \frac{∂}{∂μ}ψ(u))^{t}. El límite en (14) es equivalente a

\frac{\sqrt{T} (ψ̂(u) − ψ(u))}{σ_{R}(θ, u)} ⟶^{d} N_{1}(0, 1), T → ∞,

donde

σ_{R}(θ, u) = [ψ′(u)^{t}I_{1}^{−1}(θ)ψ′(u)]^{1/2}.

De lo anterior y de (9), se sigue

σ_{R}(θ, u) = \frac{μ}{c}e^{−(1/μ−λ/c)u}\sqrt{λ (1 + \frac{λu}{c})^{2} + λ (1 + \frac{u}{μ})^{2}}.  (15)

Por lo tanto, un intervalo asintótico del 100(1 − α)% de confianza para la probabilidad de ruina está dado por

ψ̂(u) ± z_{1−α/2} \frac{σ_{R}(θ̂, u)}{\sqrt{T}}.  (16)

Rev.Mate.Teor.Aplic. (ISSN print: 1409-2433; online: 2215-3373) Vol. 29(2): 239–260, Jul–Dec 2022

Observe que el límite superior (LS) del intervalo anterior se puede escribir como

LS = z_{1−α/2} \frac{K(θ̂, u)}{\sqrt{T}}e^{−r^{∗}u},

donde por la Proposición 3.1, K(θ̂, u) satisface

\frac{K(θ̂, u)}{\sqrt{T}} ⟶^{c.s.} 0.

De aquí, con probabilidad uno, para 0 < δ < 1 y T suficientemente grande, LS < δe^{−r^{∗}u} < e^{−r^{∗}u}. En otras palabras, para T suficientemente grande, con un nivel aproximado del 100(1 − α)% de confianza la probabilidad de ruina se puede acotar con una cantidad menor que la cota superior de Lundberg.

En general, debido a que no se conoce de manera explícita la probabilidad de ruina ψ(u) y por lo tanto no se conocen las derivadas respecto a los parámetros que la definen, no es posible calcular σ_{R}(θ̂, u). En [9] se propone una aproximación para ψ′(u) y una forma de obtener intervalos aproximados de confianza. Bajo condiciones generales, el Teorema 7 en [9] establece lo siguiente

\frac{\sqrt{T} (ψ̂(u_{T}) − ψ(u_{T}))}{σ^{∗}(θ̂, u_{T})} ⟶^{d} N_{1}(0, 1), T → ∞,

donde {u_{T} }_{T ≥0} es una sucesión de números reales tal que, cuando T → ∞,

u_{T} → ∞, \frac{u_{T}}{\sqrt{T}} → 0 y

σ^{∗}(θ, u) := [D_{θ}^{t}I_{1}^{−1}(θ)D_{θ}]^{1/2}ue^{−r^{∗}u},

con I_{1}^{−1}(θ) = T I_{T}^{−1}(θ), D_{θ} = \frac{C_{θ}H_{1}}{E[H_{0}]}, donde la constante C_{θ}, H_{1} satisfacen

C_{θ} = \frac{c − λμ}{λM′_{Y} (r^{∗}) − c}, H_{1} = \frac{1}{c} ∫_{0}^{∞} e^{r^{∗}z} \frac{∂}{∂θ} [ λ(1 − F (z))] dz

y H_{0} es una variable aleatoria con función de distribución dada por

F_{H_{0}}(x) = \frac{1}{c} ∫_{0}^{x} e^{r^{∗}z}λ(1 − F (z))dz.

De esta manera, el intervalo asintótico de 100(1 − α)% de confianza obtenido para la probabilidad de ruina en [9] está dado por

ψ̂(u_{T}) ± z_{1−α/2} \frac{σ^{∗}(θ̂, u_{T})}{\sqrt{T}}.  (17)

En el caso de las reclamaciones exponencialmente distribuidas, es posible verificar que se satisfacen las condiciones dadas mencionadas en el Teorema 7 en [9] y con algunos cálculos adicionales (ver Apéndice), se obtiene

C_{θ} = \frac{λμ}{c}, H_{1} = (\frac{1}{λ}, \frac{c}{λμ^{2}} )^{t} , F_{H_{0}}(u) = 1 − e^{−(λ/c)u}, D_{θ} = (\frac{λμ}{c^{2}}, \frac{λ}{cμ})^{t} .

De lo anterior y de (9), se verifica

σ^{∗}(θ, u) = \frac{ue^{−r^{∗}u}}{c} \sqrt{λ (\frac{λ^{2}μ^{2}}{c^{2}} + 1)}.

Observación 3.6 Se satisface σ^{∗}(θ, u) ∼ σ_{R}(θ, u), cuando u → ∞. En efecto,

(\frac{σ^{∗}(θ, u)}{σ_{R}(θ, u)})^{2} = (\frac{u}{μ})^{2} (\frac{λ^{2}μ^{2}}{c^{2}} + 1) [(1 + \frac{λu}{c})^{2} + (1 + \frac{u}{μ})^{2}]^{−1}

= [\frac{λ^{2}}{c^{2}} + \frac{1}{μ^{2}} ] [(\frac{1}{u} + \frac{λ}{c})^{2} + (\frac{1}{u} + \frac{1}{μ})^{2}]^{−1}

→ 1,

cuando u → ∞.

4 Ejemplo numérico

Con la intención de aplicar e ilustrar los resultados establecidos en la sección anterior, se realizó un análisis de simulación del modelo de Cramér-Lundberg con los siguientes parámetros:

Parámetro Valor
λ  100
μ  0.80
u  10
c  100

Los valores de λ, μ y c satisfacen la condición de beneficio neto porque

\frac{λμ}{c} = \frac{(100)(0.80)}{100} < 1.

La simulación fue realizada con el software estadístico R. En la siguiente tabla se ilustra este proceso.

1. Generar n_{1} variables aleatorias exponenciales con la función rexp(n_{1}, λ) : T_{1}, T_{2}, ..., T_{n_{1}}.

2. Determinar el valor de N_{T} asignando el valor n_{2} que cumpla T_{1} + ... + T_{n_{2}} ≤ t < T_{1}, ..., T_{n_{2}+1}.

3. Generar n_{2} variables aleatorias exponenciales con la función rexp(n_{2}, μ^{−1}) : Y_{1}, ..., Y_{n_{2}}.

4. Calcular S_{N_{T}} tomando los valores Y_{1}, ..., Y_{n_{2}}.

Con los datos simulados se calculan los estimadores de los parámetros con las fórmulas obtenidas en (7):

Estimador Valor
N_{T}  50,381
λ̂  100.7620
μ̂  0.8061
r̂^{∗}  0.2329

En la siguiente tabla se presentan la probabilidad de ruina ψ̂(u_{T}), las desviaciones estándar estimadas σ_{R}(θ̂, u_{T}), σ^{∗}(θ̂, u_{T}) para diferentes valores de T así como el intervalo de confianza al 95% de la probabilidad de ruina obtenido mediante el método delta. Cabe mencionar que en la ecuación (15) no se requiere la sucesión {u_{T} } indicada en el Teorema 7 de [9], sin embargo, con el fin de comparar σ_{R} con σ^{∗} consideramos la sucesión \{u_{T} \} definida por u_{T} = u ∗ \sqrt[8]{T} con T ∈ [100, 500]:

Tiempo ψ(u_{T}) ψ̂(u_{T}) σ_{R}(θ̂, u_{T}) σ^{∗}(θ̂, u_{T})  IC (95%)
100 0.00938 0.01039 0.0240  0.0295 (0.0056, 0.0151)
200 0.00627 0.00701 0.0171  0.0217 (0.0046, 0.0093)
300 0.00487 0.00548 0.0138  0.0178 (0.0039, 0.0070)
500 0.00348 0.00394 0.0104  0.0136 (0.0030, 0.0048)

En la Figura 2se presenta la gráfica de las probabilidades de ruina verdaderas con los parámetros propuestos al inicio de esta sección, así como las estimaciones de las probabilidades de ruina y los intervalos de confianza obtenidos con los datos simulados.

La estimación de la probabilidad de ruina tiene una mejor aproximación al valor verdadero conforme el valor de T aumenta, consecuencia directa de que el resultado es asintótico y razón por la cual para valores pequeños de T el intervalo de confianza incluya valores negativos.

Figura 2: Estimación de la probabilidad de ruina.

Finalmente, mediante simulación Montecarlo se realizaron 10,000 simulaciones con los siguientes parámetros

Parámetro Valor
λ  100
μ  0.80
u  10
c  100
T  20,000

obteniendo una cobertura para la probabilidad de ruina igual a 0.954 y los siguientes resultados:

Parámetro Estimación Error cuadrático medio
λ̂  99.9984  0.0046
μ̂  0.7999  3.3038 x10^{−7}
ψ̂(u)  0.06565  6.2291 x10^{−7}

siendo el valor verdadero de la probabilidad de ruina: 0.0657801.

5 Conclusiones

Como fue mencionado al principio, el artículo presenta de manera detallada el procedimiento para hallar los estimadores de máxima verosimilitud de los parámetros que definen al proceso de riesgo clásico con reclamaciones exponenciales, al igual que la estimación puntual y por intervalos de la probabilidad de ruina. En lo referente a lo computacional, se calcularon las estimaciones por intervalo de la probabilidad de ruina a través del tiempo y los resultados se ilustran en la Figura 2, usando la información de la simulación realizada en la Sección 4. Los intervalos calculados con la desviación estándar aproximada (17) son ligeramente más amplias que los intervalos de confianza calculados con la desviación estándar verdadera (16). Asimismo, mediante simulación Montecarlo se realizaron 10,000 simulaciones con los parámetros λ, μ, u y c dados al inicio de la sección y considerando T = 20, 000, se obtuvo una estimación puntual prácticamente igual al valor verdadero de la probabilidad de ruina. Por otro lado, la cobertura para la probabilidad de ruina con las simulaciones realizadas fue 95.4%. De esta forma, se considera que el método delta aproxima de manera adecuada los intervalos de confianza de la probabilidad de ruina para el caso de reclamaciones individuales exponenciales.

Durante el desarrollo del ejemplo numérico, se observó que para la sucesión {u_{T} } dada, se requieren valores de T grandes para obtener mejores estimaciones de los diferentes parámetros estudiados, al igual que obtener una buena aproximación para la probabilidad de cobertura del intervalo asintótico para la probabilidad de ruina. Establecer algún criterio que permita sugerir algún valor de T a partir del cuál se tengan buenas estimaciones para las cantidades de interés, es un problema que puede ser estudiado a futuro, considerando la distribución exponencial (o cualquier otra distribución) para las reclamaciones individuales.

Apéndice

Cálculo de σ^{∗}(θ, u)

Se tiene la fórmula

σ^{∗}(θ, u) := [D_{θ}^{t}I_{1}^{−1}(θ)D_{θ}]^{1/2}ue^{−r^{∗}u},

donde I_{1}^{−1}(θ) = T I_{T}^{−1}(θ), D_{θ} = \frac{C_{θ}H_{1}}{E[H_{0}]}, H_{0} es una variable aleatoria con

función de distribución dada por

F_{H_{0}}(x) = \frac{1}{c} ∫_{0}^{x} e^{r^{∗}z}λ(1 − F (z))dz,

la constante C_{θ} y el vector H_{1} satisfacen

C_{θ} = \frac{c − λμ}{λM′_{Y} (r^{∗}) − c}, H_{1} = \frac{1}{c} ∫_{0}^{∞} e^{r^{∗}z} \frac{∂}{∂θ} [ λ(1 − F (z))] dz.

Ahora bien, si Y es una variable aleatoria exponencial con media μ, entonces

M′_{Y} (θ) = \frac{μ}{(1 − μθ)^{2}} .

De aquí,

C_{θ} = \frac{λμ}{c}.

Por otro lado, utilizando que F es la función de distribución acumulada de una variable aleatoria exponencial con media μ, el vector H_{1} resulta

H_{1} = \frac{1}{c} ∫_{0}^{∞} e^{r^{∗}z} (\frac{∂}{∂λ}(λe^{−z/μ}), \frac{∂}{∂μ}(λe^{−z/μ}))^{t} dz = (\frac{1}{λ}, \frac{c}{λμ^{2}} )^{t} ,

mientras que la distribución de H_{0} es

F_{H_{0}}(u) = \frac{λ}{c} ∫_{0}^{u} e^{r^{∗}z}e^{−z/μ}dz = 1 − e^{−(λ/c)u}.

De la expresión anterior se observa que la variable aleatoria H_{0} tiene una distribución exponencial con media c/λ, así

E[H_{0}] = \frac{c}{λ}.

Poniendo todas las piezas juntas, obtenemos

D_{θ} = \frac{C_{θ}H_{1}}{E[H_{0}]} = (\frac{λμ}{c^{2}}, \frac{λ}{cμ})^{t} .

Por lo tanto, con ayuda de (9), se obtiene

σ^{∗}(θ, u) = [D_{θ}^{t}I_{1}^{−1}(θ)D_{θ}]^{1/2}ue^{−r^{∗}u}

= [(\frac{λμ}{c^{2}}, \frac{λ}{cμ}) (\begin{matrix}λ & 0 \\ 0 & \frac{μ^{2}}{λ} \end{matrix} ) (\frac{λμ}{c^{2}} ; \frac{λ}{cμ})]^{1/2} ue^{−r^{∗}u}

= \frac{ue^{−r^{∗}u}}{c} \sqrt{λ (\frac{λ^{2}μ^{2}}{c^{2}} + 1)}.

Rev.Mate.Teor.Aplic. (ISSN print: 1409-2433; online: 2215-3373) Vol. 29(2): 239–260, Jul–Dec 2022

Agradecimientos

Los autores agradecen a los árbitros anónimos por expresar comentarios y sugerencias que permitieron mejorar este manuscrito.

Financiamiento

No hay fondos comprometidos para el desarrollo de este artículo, salvo los salarios percibidos por los autores en sus respectivas instituciones.

Referencias

[1] S. Asmussen, Ruin Probabilities, World Scientific, Singapore, 2010. Doi: 10.1142/2779

[2] G. Blom, Harald Cramer 1893-1985, The Annals of Statistics 15 (1987), no. 4, 1335–1350. Doi: 10.1214/aos/1176350596

[3] G. Casella, R.L. Berger, Statistical Inference, Thomson Learning, Pacific Grove CA, 2002.

[4] A.C. Davison, Statistical Models, Cambridge University Press, Cambridge, 2003. Doi: 10.1017/CBO9780511815850

[5] D.C.M. Dickson, Insurance Risk and Ruin, Cambridge University Press, Cambridge, 2016. Doi: 10.1017/9781316650776

[6] J. Grandell, Aspects of Risk Theory, Springer New York, NY, 1991. Doi: 10.1007/978-1-4613-9058-9

[7] E.F.O. Lundberg, I. Approximerad framställning av sannolikhetsfunktionen. II. Återförsäkring af kollektivrisker. (I. Approximations of the probability function. II. Reinsurance of collective risks.) Ph.D. Thesis, Uppsala: Almqvist and Wiksells, 1903.

[8] J. Marsden, A. Tromba, Vector Calculus, W.H. Freeman and Company, New York, 2012.

[9] T. Oshime, Y. Shimizu, Parametric inference for ruin probability in the classical risk model, Statistics & Probability Letters 133 (2018), 28–37. Doi: 10.1016/j.spl.2017.09.020

[10] T. Rolski, H. Schmidli, V. Schmidt, J.L. Teugels, Stochastic Processes for Insurance and Finance, Wiley, New York, 1999.

[11] M.J. Schervish, Theory of Statistics, Springer, New York, 1995. Doi: 10.1007/978-1-4612-4250-5

[12] H. Schmidli, Risk Theory, Springer, Cham, Switzerland, 2017. Doi: 10.1007/978-3-319-72005-0

[13] Y.K. Tse, Nonlife Actuarial Models: Theory, Methods and Evaluation, Cambridge University Press, Cambridge, 2009. Doi: 10.1017/CBO9780511812156

[14] A.W. Vaart, Asymptotic Statistics, Cambridge University Press, Cambridge, 1998. Doi: 10.1017/CBO9780511802256