Análisis de riesgo con el modelo FAIR-U, usando simulación de Monte Carlo.

1 Presentación:

Se presenta un modelo simple para evaluar el riesgo, basado en el modelo propuesto por FAIR-U, que en muchos países se ha convertido en el estándar de facto para este tipo de problemas. Se lo aplica al análisis de un caso concreto: el riesgo de propagación indebida de información de clientes de la corporación Favorita.

Se consideran dos variables de entrada: la magnitud de la pérdida producida por un evento; y la frecuencia anual de los eventos que la producen. Y una de salida: el riesgo de pérdida anual estimado.

Como sucede en la mayor parte de los estudios sobre riesgo, no se tienen datos de experiencias anteriores sobre los cuales basarse. Se recurre, entonces, a la opinión de expertos.

1.1 Modelos:

Con propósitos comparativos, se hizo primero un modelo básico tentativo, bajo la hipótesis intuitiva de que la relación entre las magnitudes y las frecuencias es lineal e inversa, con las frecuencias mayores de las pérdidas anuales asociadas con sus menores magnitudes y viceversa. Una justificación de esta hipótesis podría ser que los contratos de menor valor (con magnitud probable de pérdida más pequeña) son mucho más numerosos que los de valores altos, que conllevan una mayor probabilidad de pérdida. El modelo requiere como datos los valores mínimo, máximo y más probable de la pérdida anual estimada; junto con las estimaciones de las frecuencias anuales correspondientes a cada caso.

Luego se implementó un modelo más realista: el modelo FAIR-U, que ensancha el espectro de las relaciones entre magnitudes y frecuencias, al establecer entre ellas asociaciones aleatorias; inmersas, eso sí, dentro de las distribuciones de frecuencias de ambas variables.

En este segundo modelo se requiere que los expertos entreguen dos ternas de datos, según la práctica utilizada desde mediados del siglo XX en la metodología PERT (Program Evaluation Review Technique). La primera terna contiene los valores mínimo, máximo y más probable de la pérdida anual estimada; mientras la segunda, los valores análogos de la frecuencia anual de las pérdidas.

Ambos modelos son no determinísticos, es decir, no persiguen obtener un único resultado numérico como respuesta a la pregunta de estimación del valor anual del riesgo - sobre el cual podrían haberse realizado luego pruebas de sensibilidad variando los supuestos iniciales. Nuestro objetivo es, en cambio, obtener un intervalo de valores dentro del cual pueda ubicarse la respuesta con una confiabilidad determinada, sustentada en la inferencia estadística, a partir de una distribución de probabilidades asociada al fenómeno que se esté estudiando.

1.2 Simulación del problema en la computadora:

La herramienta matemática fundamental utilizada en modelos de este tipo es una forma particular de los métodos numéricos de Monte Carlo, aptos para resolver cualquier problema que tenga una interpretación probabilística. Uno de sus usos frecuentes es la simulación de problemas que tienen significativa incertidumbre en sus datos, como es el caso del estudio del riesgo financiero que nos ocupa.

La simulación de Monte Carlo se utiliza aquí para evaluar el riesgo y la incertidumbre que afectarían a diferentes opciones de decisión como consecuencia de eventos de riesgo como el incumplimiento de disposiciones sobre confidencialidad en el manejo de información delicada de clientes.

Al aplicarlo, este algoritmo computacional genera una muestra de longitud determinada de números pseudoaleatorios (al ser la computadora una máquina determinística, el resultado no puede ser estrictamente aleatorio, pero sí ceñirse satisfactoriamente a las estadísticas prescritas para una secuencia aleatoria). Además, la muestra generada puede ajustarse a una determinada distribución probabilística, que debe ser especificada al algoritmo generador.

1.2.1 Distribución de probabilidad:

Para cada estudio considerado, es crucial la selección de una distribución de probabilidades adecuada. De acuerdo con FAIR, la más adecuada para modelar el riesgo es la distribución beta_PERT, que ha sido implementada en nuestros dos modelos.

Sobre la base de una terna de datos, se calcula cada distribución de frecuencias requerida utilizando la distribución beta_PERT, que en último término es una distribución normal continua, derivada de la distribución beta.

La distribución teórica beta requiere dos parámetros especificativos, generalmente llamados \(\alpha_1\) y \(\alpha_2\). Éstos carecen de significado intuitivo, pero, experimentando un poco, se puede apreciar cómo la relación entre ellos afecta a la curva de distribución.

La beta_PERT, en cambio, no utiliza esos parámetros, sino los característicos del PERT ya indicados. Se trata de una versión especial de la distribución beta, dados los valores estimados mínimo, máximo y más probable, y se usa para modelar datos de opinión experta. En nuestra implementación, la función rpert toma estos valores como parámetros, y devuelve una muestra de la distribución beta_PERT. Un cuarto parámetro opcional, el valor \(\lambda\), representa qué tan puntiaguda es la distribución generada. \(\lambda = 4\) corresponde a esa característica en la curva normal; un valor de \(\lambda\) más alto refleja una mayor confianza en la estimación más probable entregada por el experto, y tiende a estrechar el intervalo de confianza resultante.

Para poder transformar los parámetros de PERT a los de beta existe un modelo matemático, esbozado a continuación:

La distribución beta está caracterizada por dos funciones: la de densidad \(f(x)\) y la de distribución \(F(x)\). En el intervalo \(0 \le x \le 1\) estas funciones se definen como:\(f(x) = \frac{x^{(\alpha_1-1)}(1-x)^{(\alpha_2-1)}}{B(\alpha_1,\alpha_2)}\) y \(F(x) = \frac{B_x(\alpha_1,\alpha_2)}{B(\alpha_1,\alpha_2)}\). Fuera de [0,1] las dos funciones valen cero.

La distribución PERT usa el parámetro de mayor probabilidad para generar los parámetros de forma \(\alpha_1\) y \(\alpha_2\) basándose en este modelo. La función rpert escrita en R, implementa ese cálculo, y luego genera la distribución requerida.

1.2.2 Tamaño adecuado de las muestras en el proceso de simulación:

Otra cuestión importante que vale la pena considerar es la de la longitud adecuada de cada muestra que se genere en la aplicación de un método de Monte Carlo. Esto puede determinar si se establece un procedimiento iterativo o no. Si se escogen muestras pequeñas, entonces, para acercarnos al límite garantizado por el teorema central del límite, deberemos realizar suficientes iteraciones, y establecer la forma en que todos los resultados parciales se utilicen en la aproximación al límite. Si el tamaño de muestra es lo bastante grande, podemos obtener el resultado en una sola aplicación del método.

Al respecto, hay tres cuestiones que deben tomarse en cuenta:

1.- Capacidad de memoria.
Esta consideración, actualmente, es generalmente irrelevante, aún en las microcomputadoras, por su gran capacidad de memoria principal (del orden de decenas de gigabytes).

2.- Longitud del generador de números pseudoaleatorios.
Los algoritmos de generación tienen un período al cabo del cual los números comienzan a repetirse (y, obviamente, ya dejan de ser aleatorios). La longitud escogida para una muestra, entonces, no debe exceder este valor. En el caso del R, los algoritmos implementados tienen períodos de al menos \(2^{32}\), es decir, más de cuatro mil millones de números.

3.- Algoritmo de agregación y convergencia
En el caso que nos ocupa, la distribución de frecuencias puede generarse directamente con una sola iteración, usando la curva de densidad correspondiente a un número suficiente de puntos. Si se usan varias iteraciones, entonces, para ahorrar memoria (que sería la justificación para iterar con muestras pequeñas) habría que calcular y guardar sólo los parámetros requeridos de cada iteración y al final promediarlos o aplicar otra operación adecuada sobre ellos. (La agregación sucesiva de todos los puntos de cada iteración, para calcular la envolvente y al final sus parámetros, no sería posible cuando la memoria está restringida, y de cualquier manera sería sólo una forma menos eficiente de llevar a cabo el mismo procedimiento directo).

En nuestro caso, para cada aplicación de Monte Carlo en cada uno de los dos modelos, se escogió una sola muestra, con un número \(n\) de puntos igual a un millón (cada arreglo ocupa entonces aproximadamente 7.6 Mb).

2 Modelo básico:

En este modelo en particular, se requiere que los expertos entreguen datos de los valores mínimo, máximo, y más probable de la pérdida anual estimada; JUNTO CON SUS ESTIMACIONES DE LAS FRECUENCIAS ANUALES DE OCURRENCIA CORRESPONDIENTES A CADA CASO.

2.1 Implementación y pasos de ejecución:

Se ha escrito un guión en R para:

1.- Recibir los datos de opinión de expertos:

Estos son, para la variable M (magnitud de pérdida), el mínimo, el valor más probable y el máximo.
A cada uno acompañará el valor estimado por el experto de la frecuencia anual F en que esa pérdida puede presentarse.

2.- Hallar los tres valores de la variable R (Riesgo), como el producto de las respectivas M y F:

Las unidades del riesgo son, entonces, dólares anuales. Por conveniencia en los gráficos las abscisas representan miles de dólares por año.

\[ R_i = F_i \times M_i, \forall i \in \{1,2,3\} \]

3.- Generar una simulación usando un número adecuado de puntos:

Se ha escogido un millón para simular, usando Monte Carlo, la distribución de frecuencias de la variable R, de acuerdo con el modelo beta-Pert, el más conveniente para este problema.

4.- Calcular e imprimir las estadísticas respectivas:

Éstas describen la distribución de riesgo anual obtenido mediante sus parámetros: mínimo, máximo, promedio, moda, mediana, desviación estándar, sesgo y curtosis.
Además los deciles, que nos ayudan a leer la curva de distribución de frecuencias generada y los percentiles del 5% y 95% que permiten especificar el intervalo de respuesta, con nivel de confianza del 90%. 

2.2 Ejemplo de solución:

A continuación se ha resuelto un ejemplo que parte de datos reales proporcionados por expertos sobre el caso concreto que interesa a este estudio, usando una distribución beta_PERT sesgada hacia la derecha y con \(\lambda = 4\), que corresponde a un ajuste normal:

2.2.1 Datos:

Se parte de los pares (valor, frecuencia anual correspondiente), proporcionados por un experto:

Mínimo: ($4000, 1)
Más probable ($15000, 0.3)
Máximo: ($65000, 0.1)

Se calculan los estimados de riesgo R correspondientes (en miles de dólares anuales):

Mínimo: 4
Más probable: 4.5
Máximo: 6.5

2.2.2 EJECUCIÓN Y RESULTADOS:

Se aplica la simulación con esos datos:

2.2.2.1 Impresión de parámetros estadísticos:

Estadísticas de la distribución:
rpert(1000000,2.4,10,7.5,lambda=4)
==================================

Mínimo:     4.000
Máximo:     6.448
Promedio:   4.750
Moda:       4.501
Mediana:    4.691
Desv.est:   0.433
Sesgo:     0.5776
Curtosis:  2.7797

                                     DECILES:

    0%     10%     20%     30%     40%     50%     60%     70%     80%     90%     100% 

  4.00    4.22    4.35    4.47    4.58    4.69    4.82    4.96    5.13    5.36    6.45  

Intervalo de 90% de confiabilidad:[ 4.15, 5.55]

3 Modelo FAIR-U:

3.1 Implementación y pasos de ejecución:

El algoritmo ahora tiene los siguientes pasos:

1.- Recibir los datos de opinión de expertos:

Estos consisten en dos ternas de valores:    
La primera terna contiene los valores mínimo, máximo y más probable de la pérdida anual estimada; mientras la segunda, los valores análogos de la frecuencia anual de las pérdidas. 

2.- Generar una simulación usando un número adecuado de puntos:

Luego de definir un número adecuado de puntos (1000000, según lo ya discutido), se generan en cada caso estos puntos para simular, usando Monte Carlo, las distribuciones de frecuencias de la magnitud anual de pérdida(M) y de la frecuencia anual estimada (Fr).
Para que se las pueda verificar visualmente, elabora además los respectivos histogramas.

3.- Una vez obtenidas esas dos distribuciones, se procede a utilizarlas para calcular la distribución del valor del riesgo anual:

En este método es característica la evaluación del riesgo a partir de un emparejamiento aleatorio - aunque inmerso en las respectivas formas de las distribuciones ya obtenidas - de todos los valores de magnitud y frecuencia generados.  
Este es también un proceso numérico Monte Carlo. Produce una tercera distribución de frecuencias, en este caso del riesgo estimado.           
      
Se calcula la variable R (Riesgo) como el producto de las respectivas M y Fr:
Al igual que en el modelo básico, las unidades del riesgo son dólares anuales. 

\[R_i = Fr_i \times M_i, \forall i \in [1,n]\]

y luego su distribución de frecuencias y su función de densidad, junto con los gráficos correspondientes. Por conveniencia sus abscisas representan miles de dólares por año.

4.- Calcular e imprimir las estadísticas respectivas:

Éstas describen la distribución de riesgo anual obtenido mediante sus parámetros: mínimo, máximo, promedio, moda, mediana, desviación estándar, sesgo y curtosis.
Además, los deciles, que nos ayudan a leer la curva de distribución de frecuencias generada y los percentiles del 5% y 95% que permiten especificar el intervalo de respuesta, con nivel de confianza del 90%.
Finamente, se calculan las estadísticas de esta última, utilizando la función de densidad.

3.2 Ejemplo de solución:

El mismo ejemplo anterior con datos reales, esta vez de acuerdo con lo requerido por el modelo FAIR-U:

3.2.1 Datos:

Frecuencias: mín: 0.1;  más probable: 0.3 ; máxima: 1
Magnitudes: mín: $4000 ; más probable: $15000 ; máxima: $65000
(Las magnitudes se ingresan en miles de dólares, para comodidad de uso).

3.2.2 EJECUCIÓN Y RESULTADOS:

3.2.2.1 Cálculos de las dos distribuciones, de magnitudes y de frecuencias.

lambda <- 4
# Calcula distribuciones de frecuencias y magnitudes
Fr <- rpert(1000000,0.1,1,0.3,lambda)
M <- rpert(1000000,4,65,15,lambda) 
# Dibuja sus histogramas:

3.2.2.2 Cálculo de riesgo y gráficos de histograma y función de densidad.

# Calcula Riesgos:
R <- Fr * M
# Grafica el histograma y la función de densidad
dx = density(R)
hist(R, col = "green")

plot(dx, main="Estimación de riesgo con el producto F*M") 
polygon(dx, col="blue", border="red")
grid()
perc <- quantile(R, c(0.05,0.95)); ps <- as.numeric(perc)

# Marca intervalo de 90% de confiabilidad
abline(v = perc, col = "darkgrey",lwd = 2, lty = 3)
hh <- max(as.vector(dx[[2]]));hs <- .98*c(hh,hh)
abline(h = hh)
text(ps,hs,c(sprintf('%4.2f',perc[1]),sprintf('%4.2f',perc[2])),cex = 0.6,col = 'magenta')
text(mean(ps),1.025*hh,'90%',cex=0.6, col='red')

3.2.2.3 Impresión de parámetros estadísticos:

Estadísticas de la distribución de riesgos 
==========================================

Mínimo:     0.460
Máximo:    52.172
Promedio:   8.238
Moda:       4.207
Mediana:    6.872
Desv.est:   5.498
Sesgo:     1.3917
Curtosis:  5.4733

                                     DECILES:

    0%     10%     20%     30%     40%     50%     60%     70%     80%     90%     100% 

  0.46    2.65    3.69    4.68    5.72    6.87    8.21    9.87   12.11   15.74   52.17  

Intervalo de 90% de confiabilidad:[ 2.02,19.17]

4 Conclusión:

El segundo modelo se basa más fielmente en el sugerido por FAIR-U. A diferencia del modelo básico, que empareja de antemano las magnitudes con las frecuencias que, en un comportamiento lineal del fenómeno, les habrían correspondido, el segundo calcula de antemano, independientemente, en dos vectores, valores de Frecuencias (Fr) y valores de magnitudes (M), generados mediante simulación numérica por Monte Carlo, según la distribución betaPert. Luego, calcula el riesgo R por producto de los elementos correspondientes.

De esta manera se emparejan en la multiplicación valores aleatorios de magnitudes y frecuencias, sin una regla de asociación predeterminada, que sí existía en el primer modelo. Esto produce una mayor dispersión en los valores estimados de riesgo, que obviamente se refleja en un impacto económico anual mucho más alto que en el método básico, y, al no estar sometida a una regla rígida, refleja de mejor manera el comportamiento del fenómeno en la realidad.

El resultado, dado el enfoque probabilístico de ambos métodos, no aparece como un solo valor, sino como un intervalo, que tiene un nivel de confianza dado, y ha sido calculado por inferencia estadística a partir de la distribución de probabilidades del riesgo.

Para comparar ambas curvas es bueno recordar que son curvas de densidad, es decir, el área total bajo cada curva mide UNO. Como las unidades para el eje de abscisas son las mismas (miles de dólares anuales), en el eje vertical las unidades están ajustadas para lograr esa medida total unitaria del área. Por esa razón no es significativa la comparación de las ordenadas, pero sí la de las abscisas, y la de la forma de las curvas.

Como puede verse, el intervalo respuesta del segundo método (en abscisas) es apreciablemente mayor; es decir, da resultados más amplios (está por el lado de la seguridad), según ya se había previsto: el intervalo total del modelo básico es [4.00,6.447] contra [0.44, 52.59] del FAIR-U; mientras el intervalo con 90% de nivel de confianza para el modelo básico es [4.15, 5.55] versus [ 2.01,19.13] para FAIR-U.

Así pues, con un nivel de seguridad del 95% (puesto que el 5% de la cola izquierda queda incluido en la previsión), utilizando los resultados más realísticos del modelo FAIR-U, puede preverse una pérdida anual de no más de 19130 dólares al año, si es que no se toma ninguna medida para controlar el proceso.

Obviamente, todo el análisis está sustentado en la opinión entregada por los expertos. Si ésta no es correcta, los resultados tampoco lo serán.

En cuanto a la forma de las curvas, los dos parámetros que la describen son la curtosis y el sesgo. La curtosis es una medida del peso combinado de las colas izquierda y derecha en relación con el resto de la distribución (estrictamente no está relacionada con lo puntiagudo de una curva, como muchas veces se afirma). En nuestro caso, en la distribución del riesgo para el modelo FAIR-U tenemos curtosis alta (mayor que el valor 3 de la curva normal); es, entonces, una distribución “leptocúrtica”, e indica que las colas tienen más peso relativo que en la distribución normal. En cambio, para el modelo básico, la curtosis menor que 3, pero cercana a ese valor, muestra que las colas pesan relativamente un poco menos que en la curva normal (distribución “mesocúrtica”).

El sesgo es una medida de la asimetría de la distribución de probabilidad con relación a su promedio. También se lo puede definir refiriéndolo a la forma acampanada de la curva normal: un valor de sesgo positivo indica que la cúspide de la curva (moda) se desplaza hacia la izquierda con respecto a la de la campana normal (o también que la cola derecha es más larga que la izquierda); mientras que un sesgo negativo indica desplazamiento de la moda hacia la derecha (cola izquierda más larga). En nuestro estudio ambos modelos tienen distribuciones de probabilidad sesgadas positivamente, pero la asimetría es bastante más marcada en el modelo FAIR-U que en el básico.