Plantilla para modelación de dinámica de sistemas

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

Variables de Estado (Stock variables) kids youngsters adults retirees deaths crimes

Variables de flujo (flow variables) criminal acts per criminal kid per year criminal acts per criminal youngster per year criminal acts per criminal adult per year criminal acts per criminal retiree per year

Variables auxiliares endógenas (endogenous auxiliary variables) births percentage of kids with criminal behavior percentage of youngsters with criminal behavior percentage of adukts with criminal behavior percentage of retirees with criminal behavior

Parámetros de simulación (variables en la frontera del sistema o exogenous auxiliary variables) annual fertility rate of adults annual fertility rate of youngsters criminal acts by others per year

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

Chunk con todos los elementos para modelar el caso.

library(deSolve)
## Warning: package 'deSolve' was built under R version 4.2.3
criminal.system <- function(t, state, parameters) {
  with(as.list(c(state,parameters)), {
    #Endogenous auxiliary variables
    births <- adults*annual.fertility.adults + youngsters*annual.fertility.youngsters # births per year
    criminal.kids <- kids * 5/100
    criminal.youngsters <- kids * 50/100
    criminal.adults <- kids * 60/100
    criminal.retirees <- kids * 10/100
      
    #Flow variables
    crimActs.kid <- 2 * criminal.kids
    crimActs.adult <- 4 * criminal.youngsters
    crimActs.youngster <- 12 * criminal.adults
    crimActs.retiree <- 4 * criminal.retirees
    crimActs.other <- 6000000
      
    #State (stock) variables (d)
    dkids <- births
    dyoungsters <- kids/average.time.kid
    dadults <- youngsters/average.time.youngster
    dretirees <- adults/average.time.adult
    ddeaths <- retirees/average.time.retiree
    dcrimes <- crimActs.kid*kids + crimActs.youngster*youngsters + crimActs.adult*adults + crimActs.retiree*retirees + crimActs.other 
    
    list(c(dkids,
           dyoungsters,
           dadults,
           dretirees,
           ddeaths,
           dcrimes))
  })
}

parameters<-c( 
  annual.fertility.adults = 3/100,
  annual.fertility.youngsters = 0.3/100,
  average.time.kid = 12,
  average.time.youngster = 12,
  average.time.adult = 40,
  average.time.retiree = 15
  )

InitialConditions <- c(
  kids = 1000000,
  youngsters = 1000000,
  adults = 3000000,
  retirees = 75000,
  deaths = 0,
  crimes = 0
  )

times <- seq(0,  #initial time
             50, #end time
             1 ) #time step


intg.method<-c("rk4")

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

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