Diagramas Causales, Estructura y Comportamiento de Sistemas Dinámicos

Instrucciones y recomendaciones

  1. Cada estudiante debe entregar de manera individual su trabajo y es responsable de manera individual por la calidad de los productos entregados.
  2. La colaboración grupal es deseable durante la resolución de esta tarea, pero cada estudiante debe demostrar capacidad crítica individual en la resolución de los ejercicios descritos en este documento.
  3. La presentación de tu trabajo es muy importante. Sé conciso y claro en tu exposición, describe adecuadamente tus gráficos y empléalos para hacer más claros y contundentes tus argumentos. No escribas argumentos o ideas incompletas, concluye tus pensamientos de manera adecuada.
  4. Si empleas tablas y/o gráficos, asegúrate de incluir notas para cada uno de ellos. Cada elemento gráfico en tu exposición debe explicarse por sí mismo. Toma como referencia las anotaciones que Sterman hace en figuras y tablas en el libro de texto.
  5. Revisa cuidadosamente cada inciso, responde a la totalidad de las preguntas.
  6. No te limites a la información presentada en cada caso. Si consideras que los casos son incompletos busca de manera individual más información y referencia adecuadamente tus fuentes.
  7. Es altamente recomendable consultar el libro de texto para resolver los ejercicios descritos en esta tarea.

Problema 1: Iniciativa de reforma constitucional en materia eléctrica (30 puntos)

Descripción del caso:

En septiembre de 2021, el Ejecutivo presentó una Iniciativa de Reforma Constitucional al Sector Eléctrico, la cual propone un sistema eléctrico nacional donde el Estado recupera su conducción, a través de la CFE, la cual se convierte en el organismo del Estado, responsable de la planeación y control del sistema, autónomo en el ejercicio de sus funciones y en su administración. También se propone que CFE preserve la seguridad energética, la autosuficiencia energética y el abastecimiento continuo de energía eléctrica a toda la población. Adicionalmente, se establece que la electricidad sea un área estratégica a cargo del Estado; incorporando la generación, conducción, transformación, distribución y abastecimiento de la energía eléctrica como procesos indivisibles.

La Iniciativa eliminaría el sistema de mercado planteado en 2013, pues todos los permisos de autoabasto, generación eléctrica y contratos de compra-venta de energía serían cancelados y los CELs quedarían abolidos. Asimismo, se propone que las plantas de CFE generen el 54% del suministro eléctrico nacional.

Revisa la nota completa en: https://ciep.mx/iniciativa-de-reforma-constitucional-en-materia-electrica-potenciales-consecuencias-en-las-finanzas-publicas/

Preguntas del caso:

1.1. Toma la perspectiva del poder legislativo para evaluar la iniciativa. Define el problema que representa esta reforma para el gobierno mexicano. Recuerda que la descripción del problema debe enfatizar la existencia de una brecha entre la situación actual del país en materia eléctrica y el estado futuro deseable. Investiga y analiza acerca de las políticas que se podrían implementar para mitigar o evitar efectos negativos de esta reforma. (5 puntos, producto: definición del problema y análisis de políticas)

La iniciativa propuesta representa una reversión significativa de las políticas energéticas implementadas en años anteriores. Mientras que la reforma energética de 2013 buscaba fomentar la competencia y la inversión privada en el sector energético, la nueva propuesta apunta a concentrar el control y la operación del sistema eléctrico en la Comisión Federal de Electricidad (CFE). Esta discrepancia entre el enfoque actual y el deseado genera incertidumbre sobre el rumbo futuro del sector eléctrico y sus implicaciones económicas y ambientales.

El futuro deseable apunta una mejor gestión de las finanzas públicas, compuesta por mercados competitivos que generen utilidades a CFE y los competidores privados responsables de proveer el suministro energético al país al mismo tiempo que se optimiza la utilización de recursos naturales y se transiciona a alternativas con menor impacto ambiental y social. No obstante, la reforma propuesta compromete el cumplimiento de dichos objetivos por todas las consecuencias legales, económicas y ambientales implicadas a esta decisión. La implementación de la iniciativa conllevaría costos significativos para el gobierno mexicano. Estos costos incluyen indemnizaciones a inversionistas privados por la cancelación de contratos y permisos, aumentos en los costos de operación y mantenimiento de las plantas de generación de la CFE, así como mayores requerimientos de inversión pública para mantener los niveles de capacidad y generación eléctrica. Además, se esperaría un aumento en los subsidios a las tarifas eléctricas, así como costos adicionales por externalidades ambientales y de salud. De manera que no se vislumbran indicadores tangibles que señalen una mejora a las condiciones actuales.

Con el fin de poder mitigar y/o evadir los efectos adversos de la reforma en cuestión se podrían implementar algunas de las siguientes políticas:

  1. Promover la competencia y la inversión privada de manera equilibrada: En lugar de eliminar por completo el sistema de mercado y la participación privada en el sector eléctrico, el gobierno podría buscar un enfoque que combine la presencia de la CFE con la competencia y la inversión privada. Esto podría lograrse mediante la revisión y mejora del marco regulatorio para garantizar la transparencia, la competencia justa y la protección de los intereses de los consumidores sin perjudicar las inversiones y proyectos en curso.
  2. Negociación y compensación justa con inversionistas privados: En lugar de cancelar unilateralmente los contratos y permisos de generación del sector privado, el gobierno podría buscar negociaciones con los inversionistas para llegar a acuerdos equitativos de compensación. Esto podría ayudar a evitar arbitrajes internacionales costosos y proteger la imagen del país como destino confiable para la inversión extranjera.
  3. Promover la eficiencia y la diversificación energética: El gobierno podría implementar políticas que fomenten la eficiencia energética y la diversificación de la matriz eléctrica hacia fuentes renovables y menos contaminantes. Esto no solo ayudaría a reducir los costos de operación y mitigar los impactos ambientales, sino que también contribuiría a la seguridad energética a largo plazo.
  4. Fortalecer la capacidad regulatoria y de supervisión: Es fundamental fortalecer los organismos reguladores y de supervisión del sector eléctrico para garantizar un funcionamiento eficiente y transparente del mercado. Esto incluye garantizar la independencia y la capacidad técnica de los reguladores, así como la implementación efectiva de mecanismos de supervisión y cumplimiento de las normativas.

1.2. Realiza un inventario de las variables del sistema y clasifícalas según corresponda como variables de política, de desempeño o exógenas. (5 puntos, producto: inventario de variables del sistema)

Variables de Política: Reglamentación Ambiental Finanzas Públicas Incentivos para la IDE Inversión Pública Requerida Capacidad Instalada

Variables de Desempeño: Producción de Electricidad de CFE Costos de producción Competitividad del Mercado Subsidios

Variables Exógenas: Población Demanda de Electricidad Cancelación de Contratos Indemnización

1.3. Realiza un diagrama de sistema (5 puntos, producto: diagrama de sistema)

1.4. Desarrolla un diagrama causal. Agrega los procesos de retroalimentación que consideres más relevantes (10 puntos, producto: diagrama causal)

1.5. Usa la descripción de caso y tu propio conocimiento para graficar el comportamiento de por lo menos tres de las variables que has identificado (5 puntos, producto: comportamiento de variables)

knitr::include_graphics("C:\\Users\\Javier Cazares\\OneDrive\\Escritorio\\EGAP\\Materias\\Modelación de Sistemas\\Daño Ambiental.png")

knitr::include_graphics("C:\\Users\\Javier Cazares\\OneDrive\\Escritorio\\EGAP\\Materias\\Modelación de Sistemas\\Costos de Operación.png")

knitr::include_graphics("C:\\Users\\Javier Cazares\\OneDrive\\Escritorio\\EGAP\\Materias\\Modelación de Sistemas\\Competitividad del Mercado.png")

Problema 2: Creación de COLs (Collaborative Online Learning) y MOOCs (Massive Online Open Courses) (10 puntos)

Descripción del caso:

The board of directors of a university would like to know whether it is a good idea to start a “traditional” online distance learning programme in addition to the regular on-site programme. The members of the board would like to gain an understanding of the likely dynamics of the number of online students and on-site students over time. Similar programmes elsewhere show that:

Students enroll when entering the programme, and unenroll when finishing the programme or when quitting. The annual number of students that passes the programme is equal to a percentage of the total number of online students, which is lower than the percentage of on-site students passing. Students unenroll due to dissatisfaction with the programme. Dissatisfaction is mainly caused by the perceived quality of the teaching. Enrollment is directly proportional to the tuition fee, the perceived quality of the programme, and is highly dependent on word of mouth advertisement by current and former students. The more students there are in the online programme, the lower the costs per student.

The quality of the teaching/professors depends on experience, expertise, and motivation. Motivation is inversely proportional to their workload and the number of complaints. Complaints are the direct result of the deterioration of the quality. Improving the educational quality raises the cost per online student. And long term teaching overload results in loss of expertise.

A traditional online programme may attract students that would otherwise enroll in the onsite programme, and hence, cannibalize the on-site programme. The work related to the online programme would have to be done by the professors that teach the on-site programme, in addition to their high workload for the on-site programme.

But what if (i) new online courses would not cost professors any additional time after the materials have been developed, (ii) there would be no limit to the number of students that could follow the online course, (iii) there would almost be no costs per student, (iv) these cases would generate a lot of attention, and they would result in attracting brighter students in more advanced and fun courses to teach? Those are the promises of Massive Online Learning Courses or MOOCs.

Preguntas del caso:

2.1. Desarrolla el diagrama causal del caso (4 punto, producto: diagrama causal)

knitr::include_graphics("C:\\Users\\Javier Cazares\\OneDrive\\Escritorio\\EGAP\\Materias\\Modelación de Sistemas\\Diagrama - Ejercicio 2 .png")

2.2. Formula la hipótesis dinámica del comportamiento del sistema. Enfócate en el comportamiento de las variables que consideres son más relevantes. Describe gráficamente tu hipótesis dinámica (4 puntos, productos: hipótesis dinámica y gráfico que la respalde)

La hipótesis dinámica de este caso sugiere que a medida que la calidad de la enseñanza aumenta, la satisfacción del estudiante mejora, lo que a su vez aumenta la tasa de inscripción. Un mayor número de inscripciones lleva a una mayor demanda de cursos en línea, lo que puede resultar en una disminución de la calidad de la enseñanza debido a la sobrecarga de trabajo de los profesores. Esto puede provocar una disminución en la calidad percibida del programa y una disminución en la satisfacción del estudiante, lo que finalmente puede reducir la tasa de inscripción.

2.3. ¿Qué aconsejarías hacer a la junta directiva?, ¿deben invertir en estos cursos?, justifica tu respuesta (2 puntos, producto: justificación de la recomendación solicitada)

Recomendaría a la junta directiva invertir en cursos en línea, pero con precaución y considerando las implicaciones a largo plazo. Los MOOCs pueden ser una oportunidad para atraer a estudiantes más brillantes y avanzados, así como para reducir los costos por estudiante. Sin embargo, es importante abordar adecuadamente la calidad de la enseñanza y la sobrecarga de trabajo de los profesores. Se debe establecer un equilibrio entre la expansión del programa en línea y la preservación de la calidad educativa. Además, la junta directiva debe considerar cuidadosamente cómo promover y mantener la calidad del programa en línea para garantizar la satisfacción del estudiante y el éxito a largo plazo del programa. Así mismo, debe contemplar la elasticidad que pudiera llegar a tener la inscripción en dado caso de aumentar los costos educativos y por ende, decidir si apostar por un programa educativo más caro pero de mayor calidad (reteniendo con cargas adecuadas a sus profesores) o bien, mantener costos bajos a costo de baja motivación de sus profesores y en consecuencia insatisfacción y baja de algunos alumnos.

Problema 3: Combate contra los delitos de alto impacto (10 puntos)

Descripción del caso:

Suppose that burglaries are mainly committed by members of the organized crime (OC) and by occasional thieves. Burglaries by members of the OC and burglaries by occasional burglars then largely determine the total number of burglaries.

The Netherlands is such a small country and organized crime so well organized that the amount of burglaries by members of the OC depends largely on the effective chance of being caught for burglary in the Netherlands versus the chance of being caught for burglary in neighboring countries. The effective chance of being caught for burglary is proportional to the well-spent police man-hours available for burglaries and inversely proportional to the total police man-hours required for burglaries. Assume the police system is flexible but also rather inert: if total police man-hours required for burglaries increases/decreases, then so do the police man-hours available for burglaries, but with a delay.

Burglaries by occasional thieves mainly depend on the percentage of houses with opportunities for burglaries, which is seasonal (more in winter less in summer) and influenced by the relative media attention with regard to burglaries. The latter is proportional to the total number of burglaries divided by the acceptable number of burglaries.

Preguntas del caso:

3.1. Desarrolla el diagrama causal del caso (4 punto, producto: diagrama causal)

knitr::include_graphics("C:\\Users\\Javier Cazares\\OneDrive\\Escritorio\\EGAP\\Materias\\Modelación de Sistemas\\Robos.png")

3.2. Formula la hipótesis dinámica del comportamiento del sistema. Enfócate en el comportamiento de las variables que consideres son más relevantes. Describe gráficamente tu hipótesis dinámica (4 puntos, productos: hipótesis dinámica y gráfico que la respalde)

Los robos depende de los robos hechos por ladrones ocasionales y los de los miembros de la OC. En el caso de los robos por miembros de la OC, depende la probabilidad de ser atrapados que a su vez, es determinada por las horas de trabajo disponibles que tienen los policias para atender estas necesidades. Dado que las horas en que los policías pueden atender estos robos tiene un retraso esto seguirá propiciando que sucedan robos de manera exitosa. En la medida que ese retraso pueda ser mitigado esto ayudaría aumentar las probabilidades de atrapar a los ladrones que generan el robo y de la misma manera, reducir futuros robos.

3.3. Diseña una política que ayude a reducir el número de burglaries, gráfica el nuevo comportamiento del sistema al incluir tu política y comparalo con tu hipótesis dinámica sin política (2 puntos, producto: gráfico contrastando el comportamiento de la hipótesis dinámica con la política propuesta)

knitr::include_graphics("C:\\Users\\Javier Cazares\\OneDrive\\Escritorio\\EGAP\\Materias\\Modelación de Sistemas\\Robos con Politica.png")

Una medida que pudiera ayudar a mitigar que los ladrones tengan éxito en sus cometidos sería reducir el rezgo entre el suceso y la disponibilidad de los policias para atender el caso. Por lo anterior, analizar los principales horarios y zonas delictivas sería clave para ya tener un protocolo de intervención en el que se tenga elementos poilciacos más cercanos a los hechos y pueda elevarse gradualmente la tasa de ladrones atrapados posterior al robo. Así mismo, esto reducidría los incentivos de los ladrones a robar puesto que el valor esperado de robar cada vez sería menor, tendiendo a cero o negativo, de manera que prefieran no cometer su próximo robo.

Problema 4: Modelando la epidemia COVID (30 puntos)

Este problema tiene como objetivo guiarte en el desarrollo y análisis de un modelo dinámico empleando el lenguaje de programación R. El objetivo es que te familiarices con las herramientas de análisis en R y que inicies el desarrollo de tus propios modelos dinámicos. Para este efecto tomaremos como referencia el caso de la epidemia COVID-19, basándonos en el modelo básico de SARS descrito por Sterman (2013).

Descripción del caso:

Para iniciar tomaremos como referencia el diagrama “stock-flow” de Sterman (2013) mostrado en la siguiente figura.

Como primer paso haremos un listado del tipo de variables en el modelo:

Variables de Estado (Stock variables)

  • Population Susceptible to COVID
  • Population Infected with COVID

Variables de flujo (flow variables)

  • Infection Rate

Variables auxiliares endógenas (endogenous auxiliary variables)

  • Susceptible Contatcs
  • Probability of Contact with Infected Person
  • Contacts Between Infected and Uninfected People

Parámetros de simulación (variables en la frontera del sistema o exogenous auxiliary variables)

  • Contact Frequency
  • Total Population
  • Infectivity

Esta clasificación es útil para organizar el espacio de trabajo en Rstudio.

Tutorial de modelado del caso en R:

Iniciaremos por definir este modelo dinámico como una función definida por el usuario siguiendo los siguientes pasos.

#Carga la librería deSolve empleando la función library() 
library("deSolve")

Declara el espacio de trabajo como una función definida por el usuario. En este caso sólo tienes que cambiar el nombre de la función, manteniendo el template que hemos usado en otros modelos

covid.epidemic <- function(t, state, parameters) {
  with(as.list(c(state,parameters)), {
    #Endogenous auxiliary variables
      
    #Flow variables
      
    #State (stock) variables
      
    list(c())
  })
}

Ahora agregaremos las variables de estado y los flujos asociadas a ellas. Iniciaremos con la variable de estado “Population Susceptible to COVID” como se muestra en el chunk siguiente.

covid.epidemic <- function(t, state, parameters) {
  with(as.list(c(state,parameters)), {
    #Endogenous auxiliary variables
      
    #Flow variables
      
    #State (stock) variables
    dpopulation.susceptible.to.COVID<-
      
    list(c())
  })
}

Nota que en R las variables no pueden ser nombradas con espacios. De manera que esta variable de estado es nombrada como: “population.susceptible.to.COVID” nota también que las variables de estado son precedidas por el símbolo diferencial “d” para indicar que esta es una variable de estado, la sintaxis completa es: “dpopulation.susceptible.to.COVID”. Después de este paso habrás creado exitosamente la primera variable de estado del modelo.

Una práctica muy recomendable es documentar tus modelos. En R puedes hacer esto empleando el símbolo “#” y escribiendo delante de éste una breve descripción de la variable que estas representando. Además de describir la variable que estas creando es muy recomendable escribir también las unidades de medición de tu variable. El siguiente chunk muestra un ejemplo:

covid.epidemic <- function(t, state, parameters) {
  with(as.list(c(state,parameters)), {
    #Endogenous auxiliary variables
      
    #Flow variables
      
    #State (stock) variables
    #The population susceptible to COVID is equal to the 
    #population susceptible prior to the onset of the disease
    #less all of those that have contracted it
    dpopulation.susceptible.to.COVID<-(-1)*Infection.Rate #Stock units: People/time
    
    #Esta linea indica que variables se imprimen como 
    #resultado de la integración del modelo, todas las
    #variables de estado deben estar listadas
    list(c())
  })
}

Como se muestra en el stock-flow diagram, la variable “Infection.Rate” está conectada a la variable de estado como un flujo de salida (i.e. outflow variable), esta estructura se modela como se muestra en el chunk anterior. Nota que esta variable de flujo afecta con signo negativo a la variable de estado. Este signo negativo específica a esta variable de flujo como un flujo de salida.

Podemos seguir un proceso similar para modelar y documentar la segunda variable de estado “population.infected.with.COVID” como se muestra en la siguiente figura. Nota que en este caso la variable de flujo “Infection.Rate” es modelada como un flujo de entrada (i.e. inflow variable) y por esta razón no se incluye un signo negativo y nota también que la variable de estado “population.infected.with.COVID” es especificada precedida del signo diferencial “d”.

covid.epidemic <- function(t, state, parameters) {
  with(as.list(c(state,parameters)), {
    #Endogenous auxiliary variables
      
    #Flow variables
      
    #State (stock) variables
    dpopulation.susceptible.to.COVID<-(-1)*Infection.Rate #Stock units: People/time
    dpopulation.infected.with.COVID<-Infection.Rate #Stock units: People/time
    
    list(c())
  })
}

Como se muestra en el stock-flow diagram, en este modelo solo existe una sola variable de flujo que está conectada a las variables de estado. Esta variable de flujo es determinada por dos variables auxiliares: “Contacts between Infected and Uninfected People” e “Infectivity”. Para este caso especificamos la variable de flujo “Infection.Rate” como la multiplicación simple de estas dos variables auxiliares tal como se muestra en la siguiente figura. Nota que al agregar esta nueva variable también hemos incluido la documentación que describe esta variable y sus unidades de medición.

covid.epidemic <- function(t, state, parameters) {
  with(as.list(c(state,parameters)), {
    #Endogenous auxiliary variables
      
    #Flow variables
    Infection.Rate<-Infectivity*Contacts.bt.Infected.and.Uninfected.People #[people/time]
      
    #State (stock) variables
    dpopulation.susceptible.to.COVID<-(-1)*Infection.Rate #Stock units: People/time
    dpopulation.infected.with.COVID<-Infection.Rate #Stock units: People/time
    
    list(c())
  })
}

Al definir la variable de flujo “Infection.Rate” hemos agregado dos nuevas variables que debemos especificar. La variable “Contacts between Infected and Uninfected People” es una variable auxiliar endógena que es determinada por la variable “Susceptible Contacts” y la variable “Probability of Contact with Infected Person” esta interacción es especificada de la siguiente manera:

covid.epidemic <- function(t, state, parameters) {
  with(as.list(c(state,parameters)), {
    #Endogenous auxiliary variables
    Contacts.bt.Infected.and.Uninfected.People<-Susceptible.Contacts*Probability.of.Contact.with.Infected #[people/time]
      
    #Flow variables
    Infection.Rate<-Infectivity*Contacts.bt.Infected.and.Uninfected.People #[people/time]
      
    #State (stock) variables
    dpopulation.susceptible.to.COVID<-(-1)*Infection.Rate #Stock units: People/time
    dpopulation.infected.with.COVID<-Infection.Rate #Stock units: People/time
    
    list(c())
  })
}

Al definir esta nueva variable hemos agregado dos nuevas variables auxiliares endógenas que debemos definir. La primera variable “Susceptible Contacts” es determinada por la variable de estado “Population Susceptible to COVID” y por la variable “Contact Frequency” y es especificada como se muestra en la figura siguiente. Nota que en esta ocasión al emplear la variable de estado para definir otra variable endógena auxiliar no es necesario usar el símbolo diferencial antes de la variable de estado.

covid.epidemic <- function(t, state, parameters) {
  with(as.list(c(state,parameters)), {
    #Endogenous auxiliary variables
    Susceptible.Contacts<-population.susceptible.to.COVID*Contact.Frequency #[people/time]
    Contacts.bt.Infected.and.Uninfected.People<-Susceptible.Contacts*Probability.of.Contact.with.Infected #[people/time]
      
    #Flow variables
    Infection.Rate<-Infectivity*Contacts.bt.Infected.and.Uninfected.People #[people/time]
      
    #State (stock) variables
    dpopulation.susceptible.to.COVID<-(-1)*Infection.Rate #Stock units: People/time
    dpopulation.infected.with.COVID<-Infection.Rate #Stock units: People/time
    
    list(c())
  })
}

La última variable endógena auxiliar por modelar es la variable “Probability of Contact with Infected Person”. De manera similar al caso anterior, esta variable es determinada por la variable de estado “population infected with COVID” dividida por la variable exógena auxiliar “Total Population”

covid.epidemic <- function(t, state, parameters) {
  with(as.list(c(state,parameters)), {
    #Endogenous auxiliary variables
    Probability.of.Contact.with.Infected<-population.infected.with.COVID/Total.Population #dimensionless
    Susceptible.Contacts<-population.susceptible.to.COVID*Contact.Frequency #[people/time]
    Contacts.bt.Infected.and.Uninfected.People<-Susceptible.Contacts*Probability.of.Contact.with.Infected #[people/time]
      
    #Flow variables
    Infection.Rate<-Infectivity*Contacts.bt.Infected.and.Uninfected.People #[people/time]
      
    #State (stock) variables
    dpopulation.susceptible.to.COVID<-(-1)*Infection.Rate #Stock units: People/time
    dpopulation.infected.with.COVID<-Infection.Rate #Stock units: People/time
    
    list(c())
  })
}

Esto concluye la especificación de todos los elementos endógenos del modelo. El siguiente paso es especificar los valores de los parámetros (i.e. variables auxiliares exógenas), las condiciones iniciales de las variables de estado, el horizonte temporal de análisis y el método de integración para la simulación.

En nuestro modelo hemos especificado tres parámetros: “Contact Frequency”, “Total Population” e “Infectivity”.

En R podemos especificar los parámetros del modelo empleando un vector como se muestra en la siguiente figura. Nota que el nombre de los parámetros es idéntico al nombre empleando en la especificación descrita en los pasos anteriores. También nota que cada parámetro está asociado a un valor numérico único para el cual correremos el modelo de simulación. Finalmente nota que para cada parámetro se especifican sus unidades de medición.

parameters<-c(Infectivity = 0.1, # [1] dimmensionless
              Contact.Frequency = 2, # people/day
              Total.Population = 350 ) #people)

De igual manera podemos emplear un vector en R para definir las condiciones iniciales de cada variable de estado, esto se muestra en la siguiente figura. Nuevamente el nombre de las variables de estado es idéntico al empleado en la especificación de la modelo y cada variable de estado es inicializada a un valor único.

InitialConditions <- c(population.susceptible.to.COVID = 349 ,
                       population.infected.with.COVID = 1)

Ahora es necesario especificar el vector de tiempo que será usado para simular el modelo. Esta especificación se lleva a cabo como se muestra en la siguiente figura. Nota que para especificar este vector de tiempo empleamos la función seq(). Esta función crea una secuencia numérica y requiere tres parámetros: valor inicial, valor final e intervalo de crecimiento. En nuestro modelo dinámico estos tres parámetros representan el tiempo inicial de la simulación, el tiempo final de la simulación y la resolución temporal de la simulación. Nota que para cada parámetro hemos indicado la unidad de tiempo correspondiente.

times <- seq(0 , #initial time, days
             120 , #end time, days
             0.25 ) #time step, days

Aún no discutimos con detalle las propiedades de los diferentes métodos de integración que podemos emplear para simular nuestro modelo. Por lo pronto, para este ejercicio, elegiremos el método Runge-Kutta de Orden 4, en clases subsecuentes estudiaremos las propiedades de otros métodos de integración. La figura siguiente muestra la forma de especificar este método de integración.

intg.method<-c("rk4")

El último paso en el proceso de especificación del modelo es elegir las variables que serán “impresas” por la simulación. Este es un paso muy importante ya que nos permite elegir las variables que deseamos analizar. La figura siguiente muestra la forma de especificar las variables a imprimir por el modelo. Nota que esto es especificado en la última línea de código de la función contiene nuestro modelo dinámico. Nota también que en este caso hemos elegido imprimir sólo las variables de estado y que ambas están precedidas por el signo diferencial “d” y son concatenadas empleando la función c().

covid.epidemic <- function(t, state, parameters) {
  with(as.list(c(state,parameters)), {
    #Endogenous auxiliary variables
    Probability.of.Contact.with.Infected<-population.infected.with.COVID/Total.Population #dimensionless
    Susceptible.Contacts<-population.susceptible.to.COVID*Contact.Frequency #[people/time]
    Contacts.bt.Infected.and.Uninfected.People<-Susceptible.Contacts*Probability.of.Contact.with.Infected #[people/time]
      
    #Flow variables
    Infection.Rate<-Infectivity*Contacts.bt.Infected.and.Uninfected.People #[people/time]
      
    #State (stock) variables
    dpopulation.susceptible.to.COVID<-(-1)*Infection.Rate #Stock units: People/time
    dpopulation.infected.with.COVID<-Infection.Rate #Stock units: People/time
    
    list(c(dpopulation.susceptible.to.COVID,
           dpopulation.infected.with.COVID))
  })
}

Esto concluye la especificación del modelo de simulación. El siguiente paso es ejecutar este código y llevar a cabo la simulación. Corre los chunks que contienen la librería deSolve, los vectores que guardamos como parameters, InitialConditions, times e intg.method, y también corre el chunk de la función covid.epidemic.

Una vez realizado esto, en la consola de R ya están cargados tanto el modelo como todos los parámetros necesarios para llevar a cabo la simulación, pero aún no hemos generado datos de la simulación. Para hacer esto, crearemos una base de datos “out” que contiene los resultados de la simulación empleando la función “ode”. Nota que la función “ode” emplea como parámetros de entrada condiciones iniciales de las variables de estado del modelo, el vector de tiempo, la función covid.epidemic que describe nuestro modelo, las variables exógenas (i.e. parámetros) y el método de integración. Todos estos elementos los hemos definido ya en los pasos anteriores.

out <- ode(y = InitialConditions,
           times = times,
           func = covid.epidemic,
           parms = parameters,
           method =intg.method )

Si has hecho esto correctamente verás un nuevo objeto llamado “out” listado en el panel superior derecho.

El paso final es analizar gráficamente los resultados para esto emplea la función “plot”

plot(out,
     col=c("blue"))

Preguntas del caso:

4.1. Incluye tu modelo en un sólo chunk de código en el que se utilice la función plot para ver el comportamiento de las variables de estado. Toma en cuenta que el modelo que envíes será revisado en términos de las sintaxis del código, pero fundamentalmente en términos de su funcionalidad. Un modelo que no funcione al ser ejecutado será penalizado notablemente (2 puntos, productos: chunk con el modelo).

covid.epidemic <- function(t, state, parameters) {
  with(as.list(c(state,parameters)), {
    #Endogenous auxiliary variables
    Probability.of.Contact.with.Infected<-population.infected.with.COVID/Total.Population #dimensionless
    Susceptible.Contacts<-population.susceptible.to.COVID*Contact.Frequency #[people/time]
    Contacts.bt.Infected.and.Uninfected.People<-Susceptible.Contacts*Probability.of.Contact.with.Infected #[people/time]
    
    
    #Flow variables
    Infection.Rate<-Infectivity*Contacts.bt.Infected.and.Uninfected.People #[people/time]
    
    #State (stock) variables
    dpopulation.susceptible.to.COVID<-(-1)*Infection.Rate #Stock units: People/time
    dpopulation.infected.with.COVID<-Infection.Rate #Stock units: People/time
    
    list(c(dpopulation.susceptible.to.COVID,
           dpopulation.infected.with.COVID))
  })
}

parameters<-c(Infectivity = 0.1, # [1] dimmensionless
              Contact.Frequency = 2, # people/day
              Total.Population = 350 ) #people)

InitialConditions <- c(population.susceptible.to.COVID = 349 ,
                       population.infected.with.COVID = 1)

times <- seq(0 , #initial time, days
             120 , #end time, days
             0.25 ) #time step, days

intg.method<-c("rk4")

out <- ode(y = InitialConditions,
           times = times,
           func = covid.epidemic,
           parms = parameters,
           method =intg.method )

plot(out,
     col=c("blue"))

4.2. ¿Qué sucede cuando inicializas la variable de estado “Population Infected with COVID” en cero? Explica brevemente que origina el comportamiento que observas, emplea la estructura del modelo para cimentar tu argumentación (2 puntos, productos: grafico de simulación y breve descripción).

covid.epidemic <- function(t, state, parameters) {
  with(as.list(c(state,parameters)), {
    #Endogenous auxiliary variables
    Probability.of.Contact.with.Infected<-population.infected.with.COVID/Total.Population #dimensionless
    Susceptible.Contacts<-population.susceptible.to.COVID*Contact.Frequency #[people/time]
    Contacts.bt.Infected.and.Uninfected.People<-Susceptible.Contacts*Probability.of.Contact.with.Infected #[people/time]
    
    
    #Flow variables
    Infection.Rate<-Infectivity*Contacts.bt.Infected.and.Uninfected.People #[people/time]
    
    #State (stock) variables
    dpopulation.susceptible.to.COVID<-(-1)*Infection.Rate #Stock units: People/time
    dpopulation.infected.with.COVID<-Infection.Rate #Stock units: People/time
    
    list(c(dpopulation.susceptible.to.COVID,
           dpopulation.infected.with.COVID))
  })
}

parameters<-c(Infectivity = 0.1, # [1] dimmensionless
              Contact.Frequency = 2, # people/day
              Total.Population = 350 ) #people)

InitialConditions <- c(population.susceptible.to.COVID = 350 ,
                       population.infected.with.COVID = 0)

times <- seq(0 , #initial time, days
             120 , #end time, days
             0.25 ) #time step, days

intg.method<-c("rk4")

out <- ode(y = InitialConditions,
           times = times,
           func = covid.epidemic,
           parms = parameters,
           method =intg.method )

plot(out,
     col=c("green"))

En este caso la población susceptible a infección y la infectada se mantienen sin cambios ya que al no haber ningun infectado originalmente, no existirá población expuesta a la enfermedad y por ende, tampoco habrá gente contagiada.

4.3. ¿Cómo cambia la dinámica de comportamiento del modelo si inicializas esta variable de estado a un valor positivo diferente de cero? (2 puntos, productos: gráfico de simulación y breve descripción).

covid.epidemic <- function(t, state, parameters) {
  with(as.list(c(state,parameters)), {
    #Endogenous auxiliary variables
    Probability.of.Contact.with.Infected<-population.infected.with.COVID/Total.Population #dimensionless
    Susceptible.Contacts<-population.susceptible.to.COVID*Contact.Frequency #[people/time]
    Contacts.bt.Infected.and.Uninfected.People<-Susceptible.Contacts*Probability.of.Contact.with.Infected #[people/time]
    
    
    #Flow variables
    Infection.Rate<-Infectivity*Contacts.bt.Infected.and.Uninfected.People #[people/time]
    
    #State (stock) variables
    dpopulation.susceptible.to.COVID<-(-1)*Infection.Rate #Stock units: People/time
    dpopulation.infected.with.COVID<-Infection.Rate #Stock units: People/time
    
    list(c(dpopulation.susceptible.to.COVID,
           dpopulation.infected.with.COVID))
  })
}

parameters<-c(Infectivity = 0.1, # [1] dimmensionless
              Contact.Frequency = 2, # people/day
              Total.Population = 350 ) #people)

InitialConditions <- c(population.susceptible.to.COVID = 348 ,
                       population.infected.with.COVID = 2)

times <- seq(0 , #initial time, days
             120 , #end time, days
             0.25 ) #time step, days

intg.method<-c("rk4")

out <- ode(y = InitialConditions,
           times = times,
           func = covid.epidemic,
           parms = parameters,
           method =intg.method )

plot(out,
     col=c("red"))

En este caso no se ve un cambio tan drástico ya que pese a que la cantidad inicial aumenta (+1) la realidad es que el resto de paramteros no favorecen a que la infección se propague de manera mucho más acelerada.

4.4. ¿Cómo cambia la dinámica del sistema si aumenta el valor del parámetro “Contact Frequency”? ¿El valor de este parámetro modifica el valor final de la variable de estado “Population Infected with COVID”? Explica porque sí o porque no haciendo referencia a la estructura del modelo y a los resultados de la simulación (4 puntos, productos: gráfico con resultados de simulación y descripción de resultados).

covid.epidemic <- function(t, state, parameters) {
  with(as.list(c(state,parameters)), {
    #Endogenous auxiliary variables
    Probability.of.Contact.with.Infected<-population.infected.with.COVID/Total.Population #dimensionless
    Susceptible.Contacts<-population.susceptible.to.COVID*Contact.Frequency #[people/time]
    Contacts.bt.Infected.and.Uninfected.People<-Susceptible.Contacts*Probability.of.Contact.with.Infected #[people/time]
    
    
    #Flow variables
    Infection.Rate<-Infectivity*Contacts.bt.Infected.and.Uninfected.People #[people/time]
    
    #State (stock) variables
    dpopulation.susceptible.to.COVID<-(-1)*Infection.Rate #Stock units: People/time
    dpopulation.infected.with.COVID<-Infection.Rate #Stock units: People/time
    
    list(c(dpopulation.susceptible.to.COVID,
           dpopulation.infected.with.COVID))
  })
}

parameters<-c(Infectivity = 0.1, # [1] dimmensionless
              Contact.Frequency = 4, # people/day
              Total.Population = 350 ) #people)

InitialConditions <- c(population.susceptible.to.COVID = 349 ,
                       population.infected.with.COVID = 1)

times <- seq(0 , #initial time, days
             120 , #end time, days
             0.25 ) #time step, days

intg.method<-c("rk4")

out <- ode(y = InitialConditions,
           times = times,
           func = covid.epidemic,
           parms = parameters,
           method =intg.method )

plot(out,
     col=c("magenta"))

Si, en este caso la propagación se ve mucho más acelerada a medida que aumenta el contacto. Tal como sucedió en la pandemia, la principal preocupación de los gobiernos era mitigar las aglomeraciones o eventos con muchas personas ya que el contacto con otras personas era el catalizador para la propagación del virus. En sí, la cantidad de enfermos no era un problema tan grave para la expansión del virus siempre y cuando estos estuvieran aislados del contacto de los demás. En este caso, pese a seguir siendo un número bajo de enfermos iniciales el hecho de tener un mayor contacto hacer que la propagación aumente de manera exponencial.

4.5. ¿Cómo cambia el comportamiento del modelo si la variable de flujo “Infection Rate” cambia? Sigue los siguientes lineamientos para dar tu respuesta: Responde a esta pregunta describiendo brevemente los cambios que identificas al cambiar el valor de esta variable. Emplea un par de gráficos de comportamiento del modelo para dar soporte a tu respuesta (6 puntos, productos: gráfico con resultados de simulación y descripción de resultados).

covid.epidemic <- function(t, state, parameters) {
  with(as.list(c(state,parameters)), {
    #Endogenous auxiliary variables
    Probability.of.Contact.with.Infected<-population.infected.with.COVID/Total.Population #dimensionless
    Susceptible.Contacts<-population.susceptible.to.COVID*Contact.Frequency #[people/time]
    Contacts.bt.Infected.and.Uninfected.People<-Susceptible.Contacts*Probability.of.Contact.with.Infected #[people/time]
    
    
    #Flow variables
    Infection.Rate<-Infectivity*Contacts.bt.Infected.and.Uninfected.People #[people/time]
    
    #State (stock) variables
    dpopulation.susceptible.to.COVID<-(-1)*Infection.Rate #Stock units: People/time
    dpopulation.infected.with.COVID<-Infection.Rate #Stock units: People/time
    
    list(c(dpopulation.susceptible.to.COVID,
           dpopulation.infected.with.COVID))
  })
}


# Ejemplo 1
parameters<-c(Infectivity = 0.9, # [1] dimmensionless # increase this value
              Contact.Frequency = 2, # people/day
              Total.Population = 350 ) #people)

InitialConditions <- c(population.susceptible.to.COVID = 349 ,
                       population.infected.with.COVID = 1)

times <- seq(0 , #initial time, days
             120 , #end time, days
             0.25 ) #time step, days

intg.method<-c("rk4")

out <- ode(y = InitialConditions,
           times = times,
           func = covid.epidemic,
           parms = parameters,
           method =intg.method )

plot(out,
     col=c("orange"))

# Ejemplo 2
parameters<-c(Infectivity = 0.2, # [1] dimmensionless # increase this value
              Contact.Frequency = 2, # people/day
              Total.Population = 350 ) #people)

InitialConditions <- c(population.susceptible.to.COVID = 349 ,
                       population.infected.with.COVID = 1)

times <- seq(0 , #initial time, days
             120 , #end time, days
             0.25 ) #time step, days

intg.method<-c("rk4")

out <- ode(y = InitialConditions,
           times = times,
           func = covid.epidemic,
           parms = parameters,
           method =intg.method )

plot(out,
     col=c("orange"))

En este caso podemos ver como la tasa de infectividad tiene un impacto exponencial en el cual, son aumentar la tasa de 0.1 a 0.2 el tiempo en el que la población infectada ascendaría a 350 se reduciría a pt´racticamente menos de un año (previamente se aproximanaba a los 40 meses). En este caso, esta variable es exógena a las variables de políticas que pueda manejar un determinado gobierno, por lo tanto, durante la pandemia fueron tan euxhaustivos los esfuerzos por minimizar las salidas sociales y eventos masivos ya que era lo que tenían dentro de sus posibilidades, a diferencia de la infectividad del virus que fue aumentando y disminuyendo conforme emergían nuevas variantes del virus y se generaba la famosa “inmunidad del rebaño.

4.6. El modelo que has desarrollado siguiendo el tutorial anterior es demasiado simple. Brevemente critica la formulación y estructura del modelo y lista las suposiciones del modelo que consideras son irrealistas. Uno o dos párrafos son más que suficientes para responder este punto (2 puntos, productos: respuesta textual).

El modelo asume como constantes a lo largo del periodo evaluada las variables y esto limita mucho la cercanía a la realidad, dado que, las tasas de contacto e infectividad por ejemplo eran bastante cambiantes a lo largo de las diferentes etapas de la pandemia y el surgimiento de nuevas variantes del virus.

El hecho de asumir que en promedio la gente estaba en contacto con solo 1 o 2 personas es bastante irreal dado que esto en todo caso podría aplicar para la minoría de la gente que tenía trabajo remoto. El resto de la población que tiene trabajos que implican operación o acudir presencialmente tendrían un contacto mínimo con 10 o 20 personas por poner un ejemplo. Además, el modelo aplana su curva de personas susceptibles a contagio ya que asume que a lo largo del tiempo es la misma población y que esta solo se puede contagiar una vez y permanecer contagiada. Por lo tanto, el modelo omite la realidad de como se ve la población infectada asumiendo la tasa de recuperación de personas enfermas.

En el punto anterior identificaste algunas suposiciones irrealistas. Las siguientes preguntas tienen como objetivo que explores que sucede cuando se expanda el modelo para atender sus limitaciones.

Hasta el momento hemos asumido que la población se mantiene infectada con el virus COVID de manera indefinida. En epidemiologia esto se conoce como el modelo SI (i.e. Susceptible-Infectious). El modelo SI es apropiado para representar enfermedades crónicas para las que no existe una cura. Sin embargo, en el caso de muchas enfermedades infecciosas, incluyendo COVID, SARS, viruela o influenza, las personas infectadas pueden recuperarse o en los casos más lamentables morir.

El siguiente diagra ma stock-flow expande la estructura del modelo base para describir el proceso de recuperación de la población infectada con COVID. Esta expansión del modelo en epidemiología es conocida como el modelo SIR (i.e. la “R” indica “Recovery”). Sigue las instrucciones siguientes para expandir el modelo del tutorial.

El diagrama stock-flow muestra que debes agregar tres variables nuevas: una nueva variable de estado “Population Recovered from COVID”, una nueva variable de flujo “Recovery Rate” (por simplicidad no distinguiremos entre los pacientes que se recuperan y aquellos que mueren) y un nuevo parámetro “Average Duraction of Infection”.

El parámetro “Average Duraction of Infection” indica el tiempo promedio (i.e. en días) que una persona permanece infectada con el virus COVID. Los epidemiólogos estiman que la fase de infección del COVID tiene una duración promedio de 7 a 21 días. Emplea tu criterio para elegir el valor de este parámetro.

Existen muchas formas de modelar la variable de flujo “Recovery Rate” pero la especificación empleada con mayor frecuencia es la siguiente: Recovery Rate=Population Infected with COVID/Average Duration of Infectivity

Para implementar exitosamente esta estructura en el modelo es necesario que especifiques que esta nueva variable de flujo afecta también a la variable de estado existente “Population Infected with COVID” de la siguiente manera (i.e. sintaxis en R):

dpopulation.infected.with.COVID<-Infection.Rate - Recovery.Rate

También es necesario agregar nueva variable de estado “Population Recovered from COVID”, esto lo puedes lograr empleando la siguiente especificación:

dpopulation.recovered.from.COVID<- Recovery.Rate

Recuerda que al agregar una nueva variable de estado es necesario que indiques en el vector de condiciones iniciales el valor inicial de esta variable y también indicar en la última línea de código de la función “covid.epidemic” que esta nueva variable de estado será impresa por la simulación:

list(c(dpopulation.susceptible.to.COVID, dpopulation.infected.with.COVID, dpopulation.recovered.from.COVID))

Emplea esta nueva versión del modelo para responder a las siguientes preguntas:

library("deSolve")
covid.epidemic <- function(t, state, parameters) {
  with(as.list(c(state,parameters)), {
    #Endogenous auxiliary variables
    Probability.of.Contact.with.Infected<-population.infected.with.COVID/Total.Population #dimensionless
    Susceptible.Contacts<-population.susceptible.to.COVID*Contact.Frequency #[people/time]
    Contacts.bt.Infected.and.Uninfected.People<-Susceptible.Contacts*Probability.of.Contact.with.Infected #[people/time]
    
    #Flow variables
    Infection.Rate<-Infectivity*Contacts.bt.Infected.and.Uninfected.People #[people/time]
    Recovery.Rate<-population.infected.with.COVID/average.duration.of.infectivity  #[people/time]
    
    #State (stock) variables
    dpopulation.susceptible.to.COVID<-(-1)*Infection.Rate #Stock units: People/time
    dpopulation.infected.with.COVID<-Infection.Rate-Recovery.Rate #Stock units: People/time
    dpopulation.recovered.from.COVID<- Recovery.Rate #Stock units: People/time
    
    list(c(dpopulation.susceptible.to.COVID,
    dpopulation.infected.with.COVID,
    dpopulation.recovered.from.COVID)) 
  })
}

#Asumiendo 7 día de enfermedad
parameters<-c(Infectivity = 0.1, # [1] dimmensionless
              Contact.Frequency = 2, # people/day
              average.duration.of.infectivity = 7, # days
              Total.Population = 350 ) #people)

InitialConditions <- c(population.susceptible.to.COVID = 349 ,
                       population.infected.with.COVID = 1,
                       population.recovered.from.COVID = 0)

times <- seq(0 , #initial time, days
             120 , #end time, days
             0.25 ) #time step, days

intg.method<-c("rk4")

out <- ode(y = InitialConditions,
           times = times,
           func = covid.epidemic,
           parms = parameters,
           method =intg.method )

plot(out,
     col=c("blue"))

#Asumiendo 14 día de enfermedad
parameters<-c(Infectivity = 0.1, # [1] dimmensionless
              Contact.Frequency = 2, # people/day
              average.duration.of.infectivity = 14, # days
              Total.Population = 350 ) #people)

InitialConditions <- c(population.susceptible.to.COVID = 349 ,
                       population.infected.with.COVID = 1,
                       population.recovered.from.COVID = 0)

times <- seq(0 , #initial time, days
             120 , #end time, days
             0.25 ) #time step, days

intg.method<-c("rk4")

out <- ode(y = InitialConditions,
           times = times,
           func = covid.epidemic,
           parms = parameters,
           method =intg.method )

plot(out,
     col=c("blue"))

4.7. ¿De qué manera cambia el comportamiento de la epidemia una vez que agregas estas nuevas variables al modelo? (10 puntos, productos: nueva versión del modelo, gráficos con comportamiento de las tres variables de estado).

La curva deja de aplanarse tanto para las persona susceptibles al contagio como las personas infectadas. Ahora existe una campana en la distribución temporal de la enfermedad y se hace más alta conforme se aumenta la duración del tiempo de recuperación. Por lo mismo, se cae más rápidamente la curva de población susceptible a la enfermedad, ya que por momentos habrá menos gente “disponible” para enfermarse al estar actualmente enferma y en proceso de recuperación.

4.8. ¿Describe gráficamente y con un breve texto el efecto en el sistema de cambios (i.e. incremento y decremento) de las siguientes variables: “contact frequency” y “infectivity”? Enfatiza en las diferencias que percibes con respecto del comportamiento del modelo base. Usa una concisa y breve descripción textual y gráfica (2 puntos, productos: descripción de comportamiento, gráficos describiendo el comportamiento del modelo).

library("deSolve")
covid.epidemic <- function(t, state, parameters) {
  with(as.list(c(state,parameters)), {
    #Endogenous auxiliary variables
    Probability.of.Contact.with.Infected<-population.infected.with.COVID/Total.Population #dimensionless
    Susceptible.Contacts<-population.susceptible.to.COVID*Contact.Frequency #[people/time]
    Contacts.bt.Infected.and.Uninfected.People<-Susceptible.Contacts*Probability.of.Contact.with.Infected #[people/time]
    
    #Flow variables
    Infection.Rate<-Infectivity*Contacts.bt.Infected.and.Uninfected.People #[people/time]
    Recovery.Rate<-population.infected.with.COVID/average.duration.of.infectivity  #[people/time]
    
    #State (stock) variables
    dpopulation.susceptible.to.COVID<-(-1)*Infection.Rate #Stock units: People/time
    dpopulation.infected.with.COVID<-Infection.Rate-Recovery.Rate #Stock units: People/time
    dpopulation.recovered.from.COVID<- Recovery.Rate #Stock units: People/time
    
    list(c(dpopulation.susceptible.to.COVID,
    dpopulation.infected.with.COVID,
    dpopulation.recovered.from.COVID)) 
  })
}

#Asumiendo 7 día de enfermedad
parameters<-c(Infectivity = 0.2, # [1] dimmensionless
              Contact.Frequency = 2, # people/day
              average.duration.of.infectivity = 7, # days
              Total.Population = 350 ) #people)

InitialConditions <- c(population.susceptible.to.COVID = 349 ,
                       population.infected.with.COVID = 1,
                       population.recovered.from.COVID = 0)

times <- seq(0 , #initial time, days
             120 , #end time, days
             0.25 ) #time step, days

intg.method<-c("rk4")

out <- ode(y = InitialConditions,
           times = times,
           func = covid.epidemic,
           parms = parameters,
           method =intg.method )

plot(out,
     col=c("blue"))

#Asumiendo 14 día de enfermedad
parameters<-c(Infectivity = 0.1, # [1] dimmensionless
              Contact.Frequency = 4, # people/day
              average.duration.of.infectivity = 14, # days
              Total.Population = 350 ) #people)

InitialConditions <- c(population.susceptible.to.COVID = 349 ,
                       population.infected.with.COVID = 1,
                       population.recovered.from.COVID = 0)

times <- seq(0 , #initial time, days
             120 , #end time, days
             0.25 ) #time step, days

intg.method<-c("rk4")

out <- ode(y = InitialConditions,
           times = times,
           func = covid.epidemic,
           parms = parameters,
           method =intg.method )

plot(out,
     col=c("blue"))

library("deSolve")
covid.epidemic <- function(t, state, parameters) {
  with(as.list(c(state,parameters)), {
    #Endogenous auxiliary variables
    Probability.of.Contact.with.Infected<-population.infected.with.COVID/Total.Population #dimensionless
    Susceptible.Contacts<-population.susceptible.to.COVID*Contact.Frequency #[people/time]
    Contacts.bt.Infected.and.Uninfected.People<-Susceptible.Contacts*Probability.of.Contact.with.Infected #[people/time]
    
    #Flow variables
    Infection.Rate<-Infectivity*Contacts.bt.Infected.and.Uninfected.People #[people/time]
    Recovery.Rate<-population.infected.with.COVID/average.duration.of.infectivity  #[people/time]
    
    #State (stock) variables
    dpopulation.susceptible.to.COVID<-(-1)*Infection.Rate #Stock units: People/time
    dpopulation.infected.with.COVID<-Infection.Rate-Recovery.Rate #Stock units: People/time
    dpopulation.recovered.from.COVID<- Recovery.Rate #Stock units: People/time
    
    list(c(dpopulation.susceptible.to.COVID,
    dpopulation.infected.with.COVID,
    dpopulation.recovered.from.COVID)) 
  })
}

#Asumiendo 7 día de enfermedad
parameters<-c(Infectivity = 0.1, # [1] dimmensionless
              Contact.Frequency = 4, # people/day
              average.duration.of.infectivity = 7, # days
              Total.Population = 350 ) #people)

InitialConditions <- c(population.susceptible.to.COVID = 349 ,
                       population.infected.with.COVID = 1,
                       population.recovered.from.COVID = 0)

times <- seq(0 , #initial time, days
             120 , #end time, days
             0.25 ) #time step, days

intg.method<-c("rk4")

out <- ode(y = InitialConditions,
           times = times,
           func = covid.epidemic,
           parms = parameters,
           method =intg.method )

plot(out,
     col=c("blue"))

#Asumiendo 14 día de enfermedad
parameters<-c(Infectivity = 0.2, # [1] dimmensionless
              Contact.Frequency = 2, # people/day
              average.duration.of.infectivity = 14, # days
              Total.Population = 350 ) #people)

InitialConditions <- c(population.susceptible.to.COVID = 349 ,
                       population.infected.with.COVID = 1,
                       population.recovered.from.COVID = 0)

times <- seq(0 , #initial time, days
             120 , #end time, days
             0.25 ) #time step, days

intg.method<-c("rk4")

out <- ode(y = InitialConditions,
           times = times,
           func = covid.epidemic,
           parms = parameters,
           method =intg.method )

plot(out,
     col=c("blue"))

La campana de la distribución de la población infectada se sega a la derecha conforme se aumenta la tasa de infección o el número de contactos, señalando que, a medida que estas variables se incrementen el pico de la población infectada se verá de manera más anticipada y para el final la posibilidad de infección será cada vez más mínima.

Problema 5: Pandillas y carreras armamentistas (20 puntos)

Descripción del caso:

Arms races are escalation processes between two (or more) nations or parties in a conflict that ‘watch each other and [both] respond to [uncertain] arming activities of their opponent with [relatively greater] arming activities of their own’ (Bossel 2007c, p36). The amount of weapons held by one party may actually deter the other party from attacking and vice versa, resulting in a situation of escalating armed peace. This exercise is based on the escalation model described in (Bossel 2007c, Z507).

Suppose there are two gangs, gang A and gang B. Initially, the arms stock of gang A amounts to 100% of the weapons needed to destroy gang B, and the arms stock of gang B amounts to 100% of the weapons needed to destroy gang A. The arms stock of gang A only in/decreases via the arming of gang A and the arms stock of gang B only in/decreases via the arming of gang B. Suppose that the arming of both gangs depends on an autonomous arming rate due to their intrinsic interest in arms and arming –the autonomous arming rate A and autonomous arming rate B respectively– and on an arming rate relative to the expected first order arming of the adversary, that is, the arming of gang B from the point of view of gang A without consideration of the arming of gang A and vice versa. This relative arming rate of gang A –in terms of the weapons needed to destroy gang B– then equals the product of the overassessment factor of gang B arming by gang A, the arms obsolescence rate of gang A, and this arms stock of gang B minus the arms obsolescence rate of gang A times the arms stock of gang A. The same applies to gang B: the relative arming rate of gang B –in terms of the weapons needed to destroy gang A– equals the product of the overassessment factor of gang A arming by gang B, the arms obsolescence rate of gang B and the arms stock of gang A minus the arms obsolescence rate of gang B times the arms stock of gang B. Assume that the autonomous arming rate of gang A and the autonomous arming rate of gang B are both equal to 5% of the weapons needed to destroy the other gang per month and that both arms obsolescence rates equal 10% of the weapons needed to destroy the other gang per month.

Preguntas del caso:

5.1. Desarrolla el diagrama causal del caso (5 puntos, producto: diagrama causal)

5.2. Construye un modelo de dinámica de sistemas de este caso de estudio, suponiendo que la banda A sobreestima el armamento de la banda B en un 10%, es decir, overassessment factor of gang A arming by gang B es 110%, y que la banda B estima correctamente el armamento de la banda A, es decir, overassessment factor of gang A arming by gang B es 100%. Simula el modelo durante un periodo de 100 meses. (10 puntos, producto: modelo en R mostrando el ‘plot’ del comportamiento de las variables de estado)

library("deSolve")
parameters<-c(autonomous.arming.rate.A=0.05,
              autonomous.arming.rate.B=0.05,
              obsolescence.rate.gang.A=0.10, #weapons/month 
              obsolescence.rate.gang.B=0.10, 
              overassessment.gangB.arming.gangA=1.10, 
              overassessment.gangA.arming.gangB=1.00
              )
InitialConditions <- c(arms.stock.gang.A=1.0, #weapons
                       arms.stock.gang.B=1.0 #weapons
                       )
times <- seq(0 , #initial time, months 
             100 , #end time, months 
             20 ) #time step, months
 intg.method<-c("rk4")
 arm.race <- function(t, state, parameters) {
     with(as.list(c(state,parameters)), {
     #Endogenous auxiliary variables
     relative.arming.gangA <-(overassessment.gangB.arming.gangA* obsolescence.rate.gang.A* arms.stock.gang.B)- (arms.stock.gang.A* obsolescence.rate.gang.A)
     relative.arming.gangB <-(overassessment.gangA.arming.gangB* obsolescence.rate.gang.B* arms.stock.gang.A)- (arms.stock.gang.B* obsolescence.rate.gang.B)
     #Flow variables
     arming.gangA <- autonomous.arming.rate.A* relative.arming.gangA
     arming.gangB <- autonomous.arming.rate.B* relative.arming.gangB
     #State (stock) variables
     darm.stock.gang.A <- arming.gangA
     darm.stock.gang.B <- arming.gangB
     list(c(darm.stock.gang.A,
            darm.stock.gang.B))
   }) }
#Simulate model - question 5.2
out.arm <- ode(y = InitialConditions,
              times = times,
              func = arm.race,
              parms = parameters,
              method =intg.method
              )
plot(out.arm, col=c("green"))

5.3. Describe textualmente tu hipótesis dinámica de este caso. ¿Qué comportamientos identficas? (2 puntos, producto: descripción de la hipótesis dinámica)

En este caso, la hipótesis dinámica señala que la carrera armamentística entre ambos bandos tiene un ciclo de retroalimentación positiva que impulsa el crecimiento armamentistico de cada una debido por diversos motivos. Primero, a que cada una de ellas tiene su crecimiento autónomo, no obstante, dado que una de ellas tiene una sobrevaloración hacia la tasa de armamentismo del bando rival, eso provoca una mayor sensación de amenaza y mayor impulso a continuar armandose para obtener una ventaja sobre sus adversarios. Así mismo, el hecho de que haya más interacción armada entre las bandas termina por agravar el impulso de responder al otro bando incrementando su armamento.

5.4. Ahora supongamos que la banda A subestima el armamento de la banda B en un 50%, es decir, que overassessment factor of gang B arming by gang A es del 100%-50% o del 50%, y que la banda B evalúa correctamente el armamento de la banda A, es decir, que overassessment factor of gang A arming by gang B es del 100%. Cambia el o los parámetros correspondientes y vuelve a simular el modelo durante un período de 100 meses. ¿Qué comportamiento muestra la simulación?, ¿por qué? (3 puntos, producto: descripción textual y gráfica del nuevo comportamiento del sistema)

library("deSolve")
parameters<-c(autonomous.arming.rate.A=0.05,
              autonomous.arming.rate.B=0.05,
              obsolescence.rate.gang.A=0.10, #weapons/month 
              obsolescence.rate.gang.B=0.10, 
              overassessment.gangB.arming.gangA=1.10, 
              overassessment.gangA.arming.gangB=1.00
              )
InitialConditions <- c(arms.stock.gang.A=1.0, #weapons
                       arms.stock.gang.B=1.0 #weapons
                       )
times <- seq(0 , #initial time, months 
             100 , #end time, months 
             20 ) #time step, months
 intg.method<-c("rk4")
 arm.race <- function(t, state, parameters) {
     with(as.list(c(state,parameters)), {
     #Endogenous auxiliary variables
     relative.arming.gangA <-(overassessment.gangB.arming.gangA* obsolescence.rate.gang.A* arms.stock.gang.B)- (arms.stock.gang.A* obsolescence.rate.gang.A)
     relative.arming.gangB <-(overassessment.gangA.arming.gangB* obsolescence.rate.gang.B* arms.stock.gang.A)- (arms.stock.gang.B* obsolescence.rate.gang.B)
     #Flow variables
     arming.gangA <- autonomous.arming.rate.A* relative.arming.gangA
     arming.gangB <- autonomous.arming.rate.B* relative.arming.gangB
     #State (stock) variables
     darm.stock.gang.A <- arming.gangA
     darm.stock.gang.B <- arming.gangB
     list(c(darm.stock.gang.A,
            darm.stock.gang.B))
   }) }
#Simulate model - question 5.4        
parameters<-c(autonomous.arming.rate.A=0.05,
            autonomous.arming.rate.B=0.05,
            obsolescence.rate.gang.A=0.10,
            obsolescence.rate.gang.B=0.10,
            overassessment.gangB.arming.gangA=0.5,
            overassessment.gangA.arming.gangB=1.0
            )

out.arm <- ode(y = InitialConditions,
              times = times,
              func = arm.race,
              parms = parameters,
              method =intg.method )
plot(out.arm, col=c("red"))

En este caso se observa un cambio drástico en el comportamiento del sistema. Contrario al caso anterior, ahora la banda A subestima el armamento de la banda B por lo que su armamento disminuye en lugar de aumentar con el paso del tiempo. Ya que, en este caso en lugar de sentir una mayor amenaza de parte de los rivales siente un menor riesgo. De tal manera que al reducir su tasa de armamentismo el stock se va reduciendo por el obstoletismo de algunas de ellas.El grupo B también disminuye su armamento sin embargo a una tasa menos acelerada ya que este grupo si estima correctamente la tasa de armamento de los rivales.