Simulación por el Método Montecarlo

Tal como se ha señalado en los textos previos el método Montecarlo es un método fundamentalmente estocástico. Es exactamente lo opuesto a un método deterministico ( o si lo prefieres que se expresa con un modelo matemático).

Generlidades

La familia de los métodos montecarlo son especialmente interesantes cuando tenemos un modelo matemático determinístico. Sin embargo pueden ser usados aún en casos en los que no tengamos el modelo matemático y esto es la principal ventaja de los métodos estocásticos frente a los determinísticos. Utilizar Montecarlo es declarar que no tenemos en realidad la más mínima idea de como se relacional las variables del problema y apenas podemos o queremos tener un marco de referencia para ver que riesgo asumimos al decidir.

En este texto analizaremos tres problemas típicos en los que se puede partir del modelo determínistico matemático y enriqueceremos el análisis con Montecarlo. La técina es utilizada por gran cantidad de paquetes de software comercial, entre ellos los más notables son Crystal Ball de Oracle y @Risk de Palizade Software.

Los casos que veremos tienen que ver con el modelo de Wilson para lote económico, análisis del VAN / TIR y Programación por Camino Crítico. Por una cuestión de practiciad repasaremos el modelo determinístico sólo en el caso del modelo de Wilson y dejamos al lector el repaso de los otros modelos para tener un punto de referencia con el la simulación por el método Montecarlo.

Flujo de Trabajo

Para trabajar en este método es necesario contar con un dataset (conjunto de datos históricos) o en su defecto con un modelo formal estocástico que nos permita generar ese dataset.

El mecanismo de trabajo práctico implica un enriquecimiento de la información que nos ofrece Montecarlo respecto a métodos como los que nos aporta un modelo formal (determinístico).

A efectos de poder comprender mejor esto que se expresa iniciaremos con el estudio de un modelo determinístico, para luego pasar a tomar ese modelo y someterlo al análisis ampliado del Método Montecarlo.

Seguiremos la secuencia o flujo de trabajo que se detalla:

A- Modelar un diagrama de influencia

B- Obtener Datos Históticos de las Variables

C- Si no es posible armar el dataset buscar en la bibliografía modelos que contemplen las variables

D- Si no es poneible B y/o C Modelar el problema

E- Validar el modelo (alimnetar con datos históricos para obtener valores del presente)

F- Alimentar el Modelo con datos del presente y predecir el futuro

G- Simular con Montecarlo (alimentar con datos del presente alterados por ruido estadístico) y obtener un resultado estadístico

H- Determinar el pronóstico, el nivel de confianza y establecer el riesgo del negocio o problema simulado con Montecarlo

Diagrama de Influencia

En sí mismo esta técnica ya puede considerarse un modelado del problema. Es un modelado muy simplificado, pero la mayor parte de las veces que se realiza, sobre todo si se hace en una reunión, los resultados pueden ser más que sorprendentes.

En la técnica de simulación por Dinámica de Sistemas hay un software desarrollado por el MIT que realiza este tipo de diagrama. El software se llama Ventana Simulator y puede descargarse desde:

https://vensim.com

Diagrama de Influencia

Este diagrama se construye a partir del exágono ROJO en el que se tiene a la variable que se pretende controlar u optimizar. Dentro de los cuadriláteros hay algunos que representan ratios o índice de tasas, algunos controlables por la empresa y otros independientes (tal como los impuestos). Finalmente las elipses representan variables que influyen directa o indirectamente sobre el hexágono y los colores refieren a la categorización de externas (independientes, inciertas, exógenas, o no controlables); e internas (dependientes o en las que es posible influir en parte o sobre las que tenemos certidumbre).

Modelos determinísticos

Para ofrecer un ejemplo de los pasos del plujo de trabajo abordaremos un caso de gestión de inventarios.

El Caso de la Freight Forwarder

Un caso de emprendedora exitoso en la región vitı́cola de Mednoza - Argentina está ligado a una egresada de la carrera de ingenierı́a industrial que a comienzos de la década de 1990. Ante la incertidumbre de las crisis económicas decidió, sin capital alguno, lanzarse como free lance en el negocio de mover containers con vino a los puertos de Valparaı́so y San Antonio - Chile. En los años 2000 esto se conoció como el negocio de Freight Forwarder en el que todos los activos de trabajo son subcontratado.

Su empresa actuaba como intgradora de servicios entre cada eslabon de la cadena de suministros, entre ellos transportistas, despachantes de aduana, bodegas, navieras y clientes de brokers en el extranjero.

Un factor clave de éxito ( KPI ) era el tamaño de lote óptimo a transportar en cada envı́o. Es bien conocido el modelo de Wilson que data de la década de 1920, pero que hoy resulta casi una garantı́a al fracaso si se utiliza en contextos tan turbulentos como los que nos toca vivir. A pesar de ello el modelo aún se estudia y no es descartado, pues sirve de base para ser mejorado. El modelo tiene tres variables independientes y una dependiente, a saber:

Variables Independientes

* Demanda Media Anual D (en Containers por año)

* Costo de almacenamiento unitario Ca ($ por mes)

* Costo de emisión o gestión Ce ($ por orden)

Variable Dependiente

* Tamaño de Lote Mensual EOQ

La recomendación del Método Montecarlo es salir a buscar en la bibliografía y en los papers modelos determinísticos que vinculen a estas variables. Cuanto más actuales sean los modelos mayor certeza de los resultados obtenidos tendremos.

Modelo Matemático de Wilson

\[ {EOQ} = \sqrt {(\frac {2*D * Ce} {Ca} )}\]

El caso de la Bodega Torrenti operada por el Freight Forwarder

Esta bodega tiene tres destinos principales de exportación.

  • Estados Unidos demanda 500 TEUs al año
  • Brasil demanda 53 TEUs al año
  • Cadandá demand 28 TEUs al año

  • Costo de Almacenamiento U$D 2.500

  • Costo de Gestión U$D 750

Modelado numérico en R-CRAN

Diagrama Modelo Base

D1 <-500
D2 <-52
D3 <-28

D<- D1+D2+D3

Ca <- 2500
Ce <- 750

EOQ <- sqrt(2*D*Ce/Ca)

EOQ
## [1] 18.65476

La conclusión final es que la bodega debería enviar 18 containers por mes para operar en un punto que minimice sus costos logísticos totales.

Lote Económico resuelto por Montecarlo

Este método no concibe a las variables independientes como fijas, sino que las considera variables aleatorias e introduce Ruido en el modelo matemático. En esta parte deberíaos analizar el comportamento de cada variable y asociarle una distribución (función paramétrica) de probabilidad. Por sencillez supondremos en este ejemplo que las varaibles siguen una distribución normal.
En rigor deberíamos estudiar que distribución siguen los datos con los que alimentaremos el modelo. <bR R-Cran tiene por defecto definidas las siguientes distribuciones de probabilidad.

  • Normal
  • Logística
  • Uniforme
  • Gamma
  • Lognormal
  • Weibull
  • Cauchy
  • Exponencial
  • Chi-cuadrado
  • F
  • T-Student.


Puede explorar el paquete rriskDistributions que le dará una idea más clara. Generaremos una muestra al azar de 47 valores con media 40 y desvio standard 5

#    library(rriskDistributions)
#    muestra=rnorm(47,40,5)
#    muestra
#    fit.cont(muestra)

Estocatización de Variables

  • Estados Unidos demanda promedio 500 TEUs al año con una varianza de 90

  • Cadandá demanda promedio 28 TEUs con una varianza de 6 al año

  • Brasil demanda promedio 93 TEUs al año con una varianza de 8


  • Costo de Almacenamiento U$D 2.500 Varianza 1.000

  • Costo de Gestión U$D 750 varianza 250

Cálculo del EOQ Montecarlo

#Genero tres variables aleatorias normales con 
#una cantidad de muestras m igual a tres veces #el tamaño de la media 

m1 <- 500*3
D1 <- rnorm(m1,500,90)
m2 <- 93*3
D2 <- rnorm(m2,93,8)
m3 <- 28*3
D3 <- rnorm(m3,28,6)
D <- c(D1,D2,D3)
# m_d muestras de demanda
m_d <- length(D)

#Generamos las variables aleatorias de costos igual a la cantidad de muestras de demanda que tenemos.

Ca <- rnorm(m_d,2500,1000)

Ce <- rnorm (m_d,750,250)

# Preparo área de graficos 3 filas 1 columna
#par(mfcol=c(1,3))

#Grafico 1
plot(density(D), 
xlab="Demanda TEU", #Change the x-axis label
ylab="Densidad", #y-axis label
main="Muestras de Demanda 5 años")#Main title

#Grafico 2
plot(density(Ca),
xlab="Costo Almacenamiento", #x-axis label
ylab="% eventos registrados", # y label
main="Población Costos Almacenamiento")     

#Grafico 3
plot(density(Ce),
xlab="Costo Gestión Inventario", #x-axis label
ylab="% eventos registrados", # y label
main="Población Costos Emisión")     

Cálculo de muestras de lote Económico EOQ

En este caso realizaremos las mismas operaciones que antes, pero las variables independientes ahora serán matrices y el resultado será un conjunto de valores con su propia distribución de probabilidades.

EOQ <- sqrt(2*D*Ce/Ca)
## Warning in sqrt(2 * D * Ce/Ca): Se han producido NaNs
#Media del Lote Económico
mean(EOQ, na.rm=TRUE)
## [1] 16.24418
#Varianza del Lote Económico
var(EOQ, na.rm=TRUE)
## [1] 78.75197
#Desvio Estandard del Lote Económico
sd(EOQ, na.rm=TRUE)
## [1] 8.874231
plot(density(EOQ,na.rm=TRUE))

Modelo Montecarlo

Hallaremos los valores máximos y mínimos de las variables dependientes e independientes

#Índice Mínimos
which.min(Ca)
## [1] 146
which.min(Ce)
## [1] 1435
which.min(EOQ)
## [1] 1418

Valores Mínimos

Ca[which.min(Ca)]
## [1] -1007.138
Ce[which.min(Ce)]
## [1] -28.43039
EOQ[which.min(EOQ)]
## [1] 1.603287
#Máximos
which.max(Ca)
## [1] 248
which.max(Ce)
## [1] 674
#Atención excuiremos los NA
which.max(EOQ)
## [1] 1272

Riesgo del Modelo

Generación de Dataset de Inventario Simulado

Se ha visto como a partir del modelo numérico pudimos hacer el trabajo de simulación. Si tuviésemos datos históricos con solamente calcular los parámetros \(\mu\) y \(\sigma\) podríamos reporducir el experimento y la simuación En algunas oportunidades las estadísticas como las de kaggle, CEPAL, Banco Mundial, etc. nos entregan los parámetros para construir el historial. Veremos ahora cómo generar el dataset a partir de los parámetros.

Para poder realizar este tipo de simulación es necesario que tengamos un historial de los costos de gestión de un inventario. En este ejercicio estamos pensando en un historial de datos que varios años. Por caso plantemos un ejemplo que abarca desde 1950 a 2021.

Si el caso es utilizar datos reales, todos los métodos de simulación requieren que se haga un trabajo de tratamiento y depuración de la información. Tendremos que realizar cosas y tomar decisiones que implican:

    1. Remover datos de circunstancias excepcionales o capturados erroneamente.
    1. Revisar los NA
    1. Decidir como se completaran los NA (¿Serán reemplazados por el promedio anual?, ¿Se interpolará?, ¿Se tomará una aproximación basada en el mismo dato de la misma semana del años pasado incrementado por la tendencia?, etc)
    1. Reralizar una actualización a precios constantes si corresponde
    1. Realizar entrevistas y reuniones con expertos para validar
    1. Comparar con datos de cámaras o estaddísticas nacionales Inicia en 1950, finaliza en Dic 2021
    1. Repasar la consistencia metódica de los datos

Demanda

Constituiremos el primer dataset que representa la demanda de containers adquiridos por la empresa

D1 <- round(rnorm(426,300,30))
D2 <- round(rnorm(426,200,50))
D <- c(D1,D2)
hist(D)

Costo de Emisión

Ce1 <- rnorm (426,30,4)
Ce2 <- rnorm (426,25,4)
Ce <- c(Ce1,Ce2)
plot(density(Ce))

Costo de Almacenamiento

Ca1 <- rnorm(426,72,12)
Ca2 <- rnorm(426,95,8)
Ca <- c(Ca1,Ca2)
plot(density(Ca))

Recuperación de datos desde la Web

Inventario Simulado

Si no deseas generar el datasets, puedes cargarlo con load(“Wilson.Rda”)

https://themys.sid.uncu.edu.ar/ /rpalma/R-cran/Inventario_Simulado/Wilson.Rda

load ("Wilson.Rda")
library(psych)
pairs.panels(Wilson)

Lote Económico

EOQ <- sqrt(2*D*Ca/Ce)
hist(EOQ)

Este método que hemos usado utilizado nos permite crear un dataset semejante al que podríamos obtener si capturamos los datos de la operación de un negocio formal.

Tengamos en cuenta que todos los pasos anteriores no son importantes para el trabajo que vamos a realizar. De hecho la formula del lote económico de Wilson (EOQ) fue creado en 1912. Tiene más de 100 años y en la actualidad , como modelo, garantiza el fracaso en los cotos de gestión de inventario.

Discutiremos en clase más sobre este tema.

Guardaremos este archivo y a partir de el realizaremos una simulación por el método Montecarlo.

Grabar Archivo

Guardaremos en archivo en la carpeta que nos indica el comando getwd(). Lo haremos utilizando los formatos csv (que ya conocemos) y el formato propio de R-Cran denominado .Rda

getwd()
## [1] "/media/rpalma/OS/AAA_Datos/2022/Posgrado/R-Cran/Inventario_Simulado"
Wilson <- cbind(D,Ca,Ce,EOQ)
save (Wilson, file="Wilson.Rda")
write.csv(Wilson,file="Wilson.csv")

Análisis exploratorio

Realizaremos un análisis exploratorio de cada componente del modelo. Nuestro modelo tiene una variable Dependiente llamada EOQ, que representa la cantidad de contenedores pedidos para satisfacer la demanda. y tres variavles independientes:

  • Costo de Almacenamiento Ca
  • Costo de sostenimiento o gestión del inventario Ce
  • Demanda anual D

Costo de almacenamiento

summary(Ca)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   41.27   70.97   86.29   83.36   95.45  123.72

Valores máximos y mínimos

which.max(Ca)
## [1] 811
which.min(Ca)
## [1] 82

¿Que valor tiene el costo mínimo?

Ca[3]
## [1] 75.85777

¿Qué valores están por arriba del tercer cuanti?

which(Ca> 90) 
##   [1]   1  17  43  45  55  59  73  89  92  98 112 138 153 168 173 195 209 229
##  [19] 236 244 282 292 321 325 329 333 334 350 358 359 362 396 397 427 429 430
##  [37] 431 432 433 434 436 437 439 440 441 443 444 445 446 447 448 449 450 451
##  [55] 452 456 458 459 460 461 462 463 464 465 467 468 470 472 474 475 477 479
##  [73] 480 481 483 484 486 488 489 491 492 494 496 497 498 499 500 501 503 506
##  [91] 509 510 511 512 514 515 516 517 518 520 521 523 524 525 526 527 528 529
## [109] 530 531 532 533 534 535 536 537 538 539 540 542 543 544 545 546 547 548
## [127] 549 551 552 554 555 556 557 558 560 561 562 563 565 569 570 571 572 573
## [145] 574 575 579 581 582 584 585 588 589 590 591 592 593 594 596 599 600 601
## [163] 602 603 605 607 608 609 610 611 614 616 617 618 621 623 624 625 626 627
## [181] 628 629 630 631 632 633 634 635 636 637 638 639 640 641 642 643 644 646
## [199] 647 648 649 650 651 652 653 656 657 658 659 660 661 662 663 666 668 669
## [217] 671 672 673 674 676 677 678 681 682 683 684 685 686 687 688 689 690 691
## [235] 692 693 694 695 696 698 699 700 701 703 704 705 706 707 709 710 712 713
## [253] 715 716 717 718 719 721 722 723 724 726 727 731 733 735 737 738 740 743
## [271] 745 747 748 749 750 752 753 754 755 756 757 758 760 761 763 764 766 767
## [289] 769 770 772 773 774 775 776 777 779 780 781 782 784 785 786 787 788 789
## [307] 790 791 792 793 794 795 796 797 799 800 801 802 803 804 805 807 808 810
## [325] 811 812 813 814 815 816 817 818 819 820 821 822 823 824 825 828 830 832
## [343] 833 834 835 836 837 838 839 840 842 844 845 846 847 849 850 851 852

Histogramas

Evolución de costos de almacenamiento

plot(sort(Ca), main ="costos ordenados", type="l")

DataExplorer

DataExplorer simplifica y automatiza el proceso EDA y la generación de informes. El paquete escanea automáticamente cada variable que realiza el perfilado de datos y ofrece varias funciones útiles para generar diferentes gráficos en características discretas y continuas.

Es importante destacar que el paquete permite generar un informe HTML completo sin esfuerzo al invocar la útil función create_report en un conjunto de datos (por ejemplo, create_report(df)). Es posible pasar argumentos adicionales, como una variable de respuesta y, para agregar varios análisis bivariados al informe:

Busca en la carpeta en la que has guardado este archivo. La biblioteca habrá creado un archivo llamado Reporte.html Las líneas están comentadas por que este comando no se puede correr dentro de un archivo Marckdown

# library(DataExplorer)
# Wilson %>% create_report(output_file = "Report")