Simulation - system dynamics (R).ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Simulations: system dynamics\n",
    "\n",
    "System dynamics models work on stocks and flows between the stocks.\n",
    "To get an idea about this, see how [MONIAC](https://www.youtube.com/watch?v=rVOhYROKeu4) works.\n",
    "See the documentation of [deSolve](https://cran.r-project.org/web/packages/deSolve/index.html)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "library(deSolve)\n",
    "library(tidyverse)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## time starts from 0 and continues until 20\n",
    "\n",
    "time <- seq(0, 20, by=1)\n",
    "\n",
    "## we have three stocks (people healthy, people sick and people killed)\n",
    "\n",
    "stocks <- c( healthy = 900, sick = 100, killed = 0 )\n",
    "\n",
    "## model has parameters, such as infection rate amd kill rate\n",
    "\n",
    "auxs <- c( infection = .05, dead = 0.05 )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## here we define how each step is executed\n",
    "\n",
    "model <- function(time, stocks, auxs){\n",
    "  with(as.list(c(stocks, auxs)), {\n",
    "      \n",
    "    # we can modify parameters over time\n",
    "    # if( time >= 10 ) {\n",
    "    #     infection <- .00\n",
    "    # }\n",
    "\n",
    "    new_people_sick <- sick * infection\n",
    "    new_dead_people <- sick * dead\n",
    "      \n",
    "    healthy_change <- - new_people_sick\n",
    "    dead_change <- new_dead_people\n",
    "    sick_change <- - new_dead_people + new_people_sick\n",
    "\n",
    "    ## we update all three stocks and parameters\n",
    "    return (\n",
    "        list(\n",
    "            c( healthy_change, sick_change, dead_change ),\n",
    "            auxs\n",
    "        )\n",
    "    )\n",
    "  }) \n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## this runs the model\n",
    "data <- data.frame( dede( y = stocks, times = time, func = model, parms = auxs, method = \"lsodar\"))\n",
    "head( data )\n",
    "tail( data )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df <- data %>% \n",
    "  pivot_longer(-time, names_to = 'variable')\n",
    "\n",
    "\n",
    "ggplot(df, aes(time, value, color = variable))+\n",
    "  geom_line(size =1.25)+\n",
    "  theme_minimal()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercises\n",
    "\n",
    "1. Based on the model equation, draw a picture about the model (similar to Figure 6.2)\n",
    "1. Change the infection rate to 1.10 and dead rate to 0.05 and plot the results.\n",
    "1. Increase the time span by 10 more days.\n",
    "1. Examine how a stronger intervention (decreasing the infection rate more) impacts to number of patients and deadths. What do you see?\n",
    "1. Change the time of the intervention, what changes can you detect?\n",
    "1. Examine the timespan to 300 days. Why does the model break? Fix it so that it does not!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "R",
   "language": "R",
   "name": "ir"
  },
  "language_info": {
   "codemirror_mode": "r",
   "file_extension": ".r",
   "mimetype": "text/x-r-source",
   "name": "R",
   "pygments_lexer": "r",
   "version": "4.2.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}