Tarea 3 - Modelación de sistemas

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: Escasez de minerales y metales (30 puntos)

Problema 14.15 del libro Small System Dynamics Models for Big Issues: Triple Jump towards Real-World Complexity de Erik Pruyt.

Descripción del caso:

From Extraction and Processing to Production of Goods

Rare Earths Metals (REM) are often used in very small quantities in modern appliances or applications. Assume that the REM in goods initially amounted to 2000000 ton in the year 2010. REM in goods are lost after an average lifetime in goods of some 15 years (recycling of REM is (currently) not feasible). The use of REM in the production of goods depends on the demand for REM or the available supply of REM –initially equal to 250000 ton– if the available supply of REM is smaller than the demand for REM. The available supply of REM only increases through processing of REM which in turn follows the real REM extraction.

Preguntas del caso, parte 1:

1.1. Para esta primera parte del caso construye el modelo de dinámica de sistemas en R y simula el modelo por un periodo de 80 años, tiempo inicial año 2000, tiempo final año 2080.

library("deSolve")

REM  <- function(t, state, parameters) {
  with(as.list(c(state,parameters)), {
    #Endogenous auxiliary variables
    processing.of.REM <- real.REM.extraction
    
    #Flow variables
    use.of.REM.in.the.production.of.goods <- max(0, demand.for.REM - available.supply.of.REM) 
    
    
    #State (stock) variables
    ddemand.for.REM <- use.of.REM.in.the.production.of.goods 
    davailable.supply.of.REM <-  processing.of.REM - use.of.REM.in.the.production.of.goods
    dREM.in.goods <- use.of.REM.in.the.production.of.goods/average.lifetime.in.goods
    
    list(c(ddemand.for.REM,
           davailable.supply.of.REM,
           dREM.in.goods))
  })
}

parameters<-c(average.lifetime.in.goods = 15,
              real.REM.extraction= 0) # years

InitialConditions <- c(REM.in.goods= 2000000, # ton in the year2010
                       available.supply.of.REM = 250000,
                       demand.for.REM = 0) #tons
                       
            

times <- seq(2000 , #initial time, years
             2080, #years, 
             1) #time step, years

intg.method<-c("rk4")

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

Demand and Supply of REM

The demand for REM –initially equal to 100000 tons in the year 2000– increases in principle by means of an increase of the demand for REM and decreases through a decrease of demand through price elasticity of demand or through substitution losses.

The increase of the demand for REM is simply the product of the economic growth rate and the demand for REM. Suppose for reasons of simplicity that the economic growth rate was equal to 3% until the year 2009, that it fell to -10% in 2009, and that it jumped to 8% in 2010. Assume that the economic growth rate remains constant at 5% from 2011 on.

The decrease of demand through price elasticity of demand could be modeled as:

(-1) * price elasticity of demand * demand for REM * ( (1/relative price - 1/relative price of the previous year) / (1/relative price of the previous year))

With the relative price equal to the product of the average REM extraction costs, the scarcity price effect and (1 + normal profit margin). Suppose that the price elasticity of demand is 10% and that the normal profit margin is 15%.

If the relative price is greater than the relative price of the cheapest substitute then the substitution losses increase to the product of the demand for REM, (1/relative price), and the difference between the relative price and the relative price of the cheapest substitute. The substitution losses due to price substitution effects are in other words non-negative. Suppose that the relative price of the cheapest substitute is constant at 100 (in other words 100 times the normal price of REM).

Suppose also that the scarcity price effect amounts to 100 if the supply demand ratio is 0, 10 if the supply demand ratio is 0.55, 1 if the supply demand ratio is 1.1, 0.75 if the supply demand ratio is 2.2, 0.5 if the supply demand ratio is 11, and 0.2 if the supply demand ratio is 22. The supply demand ratio is equal to the available supply of REM divided by the demand for REM.

Preguntas del caso, parte 2:

1.2. Para esta segunda parte del caso extiende el modelo de dinámica de sistemas en R y simula el modelo por un periodo de 80 años, tiempo inicial año 2000, tiempo final año 2080.

library("deSolve")

REM  <- function(t, state, parameters) {
  with(as.list(c(state,parameters)), {
    #Endogenous auxiliary variables
     
    processing.of.REM <- real.REM.extraction
    economic.growth.rate <- if (t < 2009) {
  growth_rate <- 0.03
} else if (t == 2010) {
  growth_rate <- -0.08
} else if (t == 2011) {
  growth_rate <- 0.05
} else if (t > 2009) {
  growth_rate <- -0.10
}
    supply.demand.ratio <- available.supply.of.REM / demand.for.REM
    scarcity.price.effect <- approx(c(0,0.55,1.1,2.2,11,22),#x
                             c(100,10,1,0.75,0.5,0.2),#y
                             xout=supply.demand.ratio)$y
    
    relative.price <- (average.REM.extraction.costs * scarcity.price.effect) *  (1 + normal.profit.margin)
      
    
    #Flow variables
    use.of.REM.in.the.production.of.goods <- max(0, demand.for.REM - available.supply.of.REM) 
    increase.of.the.demand.for.REM <- economic.growth.rate * demand.for.REM
    relative.price.of.the.previous.year <- relative.price
    decrease.of.demand.through.price.elasticity.of.demand <- - price.elasticity.of.demand * demand.for.REM * (((1/relative.price)-(1/relative.price.of.the.previous.year))/(1/relative.price.of.the.previous.year))
    substitution.losses <- max(0,ifelse(relative.price > relative.price.of.the.cheapest.substitute, demand.for.REM*(1/relative.price)*(relative.price-relative.price.of.the.cheapest.substitute), 0))
     
      
    
    #State (stock) variables
    ddemand.for.REM <- increase.of.the.demand.for.REM - decrease.of.demand.through.price.elasticity.of.demand - substitution.losses
    davailable.supply.of.REM <-  processing.of.REM - use.of.REM.in.the.production.of.goods
    dREM.in.goods <- use.of.REM.in.the.production.of.goods/average.lifetime.in.goods
    
    list(c(ddemand.for.REM,
           davailable.supply.of.REM,
           dREM.in.goods))
  })
}

parameters<-c(average.lifetime.in.goods = 15,
              real.REM.extraction= 0,
              price.elasticity.of.demand = .10,  # %
              normal.profit.margin = 0.15,
              relative.price.of.the.cheapest.substitute = 100,
              economic.growth.rate = 0.05,
              average.REM.extraction.costs=0) # 

InitialConditions <- c(REM.in.goods= 2000000, # ton in the year2010
                       available.supply.of.REM = 250000,
                       demand.for.REM = 100000) #tons
                       
            

times <- seq(2000 , #initial time, years
             2080, #years, 
             1) #time step, years

intg.method<-c("rk4")

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

Commissioning and Decommissioning of Extraction Capacity

Suppose that the mining industry is myopic and has limited foresight: the desired extraction capacity then equals the demand for REM. And the newly planned extraction capacity then equals the product of the profitability of REM extraction between 0 and 1 and the difference between the desired extraction capacity and the installed extraction capacity. Note: the value of the latter difference necessarily lies between 0 and the value of the installed extraction capacity.

Newly planned extraction capacity increases the extraction capacity under construction –initially equal to 60000 t/yr. The extraction capacity under construction decreases through commissioning of extraction capacity which delays the newly planned extraction capacity with a precise construction time of extraction capacity of exactly 8 years. The commissioning of extraction capacity initially equals the extraction capacity under construction divided by the precise construction time of extraction capacity.

The commissioning of extraction capacity leads of course to an increase of the installed extraction capacity, initially equal to 100000 t/yr. The installed extraction capacity decreases on the one hand through decommissioning of extraction capacity and on the other hand through decommissioning of unprofitable extraction capacity. Model the decommissioning of unprofitable extraction capacity as: installed extraction capacity * (-MIN(profitability of REM extraction,0)).

The decommissioning of extraction capacity then equals the installed extraction capacity divided by the average lifetime of extraction capacity minus the decommissioning of unprofitable extraction capacity. Note that the formula of the decommissioning of extraction capacity needs to be non-negative. Set the average lifetime of extraction capacity to 30 years.

The maximum REM extraction is equal to the installed extraction capacity. The real REM extraction normally equals the installed extraction capacity, unless the scarcity price effect is smaller than 1, then it equals the installed extraction capacity times the scarcity price effect.

Cumulation of the real REM extraction gives the cumulatively extracted REM –initially equal to 4000000 ton. The difference between the cumulatively extracted REM and the initial cumulatively extracted amount is needed to calculate the average REM extraction costs. To do so, use a function with the difference between the cumulatively extracted REM and the initial cumulatively extracted amount as argument that connects following couples: (0, 1), (2.000.000, 2), (4.000.000, 4), (6.000.000, 8), (8.000.000, 16), (10.000.000, 32), (12.000.000, 64), (14.000.000, 128), (16.000.000, 256), (18.000.000, 512).

Finally, the profitability of REM extraction is equal to the difference of the relative price and the average REM extraction costs, divided by the average REM extraction costs.

Preguntas del caso, parte 3:

1.3. Para esta tecera parte del caso extiende el modelo de dinámica de sistemas en R y simula el modelo por un periodo de 80 años, tiempo inicial año 2000, tiempo final año 2080. Modela ‘intrinsic demand’ –en otras palabras, la demanda en la ausencia de ‘decrease of demand through price elasticity of demand’ y ‘substitution losses’– y haz un indicador ‘fraction produced of intrinsic demand’ que nos permita visualizar el ‘use of REM in the production of goods’ en función de ‘intrinsic demand’ en el tiempo.

library("deSolve")

REM  <- function(t, state, parameters) {
  with(as.list(c(state,parameters)), {
    #Endogenous auxiliary variables
    
    supply.demand.ratio <- available.supply.of.REM / demand.for.REM
    scarcity.price.effect <- approx(c(0,0.55,1.1,2.2,11,22),#x
                             c(100,10,1,0.75,0.5,0.2),#y
                             xout=supply.demand.ratio)$y
    
    dif.average.REM.extraction.costs <- cumulatively.extracted.REM -  initial.cumulatively.extracted.REM
    
    average.REM.extraction.costs <- approx(c(0,2000000,4000000,6000000,8000000,10000000,12000000,14000000,16000000,18000000),#x
                             c(1,2,4,8,16,32,64,128,256,512),#y
                             xout= dif.average.REM.extraction.costs)$y
    
    processing.of.REM <- real.REM.extraction
    real.REM.extraction <- ifelse(scarcity.price.effect < 1, installed.extraction.capacity*scarcity.price.effect, installed.extraction.capacity)
    
    economic.growth.rate <- if (t < 2009) {
  growth_rate <- 0.03
} else if (t == 2010) {
  growth_rate <- -0.08
} else if (t == 2011) {
  growth_rate <- 0.05
} else if (t > 2009) {
  growth_rate <- -0.10
}

    
    relative.price <- (average.REM.extraction.costs * scarcity.price.effect) *  (1 + normal.profit.margin)
    commissioning.of.extraction.capacity <- extraction.capacity.under.construction / precise.construction.time.of.extraction.capacity
      
    
    #Flow variables
    use.of.REM.in.the.production.of.goods <- max(0, demand.for.REM - available.supply.of.REM) 
    
    increase.of.the.demand.for.REM <- economic.growth.rate * demand.for.REM
    
    relative.price.of.the.previous.year <- relative.price
    
    decrease.of.demand.through.price.elasticity.of.demand <- - price.elasticity.of.demand * demand.for.REM * (((1/relative.price)-(1/relative.price.of.the.previous.year))/(1/relative.price.of.the.previous.year))
    
    substitution.losses <- max(0,ifelse(relative.price > relative.price.of.the.cheapest.substitute, demand.for.REM*(1/relative.price)*(relative.price-relative.price.of.the.cheapest.substitute), 0))
    
    desired.extraction.capacity <- demand.for.REM
    
    profitability.of.REM.extraction <- (relative.price - average.REM.extraction.costs) / average.REM.extraction.costs
    
    newly.planned.extraction.capacity <- (min(max(profitability.of.REM.extraction, 0), 1) * min(max((desired.extraction.capacity - installed.extraction.capacity), 0), installed.extraction.capacity))
    
    decommissioning.of.unprofitable.extraction.capacity <- installed.extraction.capacity * (-min(profitability.of.REM.extraction,0))
    
    decommissioning.of.extraction.capacity <- max(0,(installed.extraction.capacity / average.lifetime.of.extraction.capacity - decommissioning.of.unprofitable.extraction.capacity))
    
    fraction.produced.of.intrinsic.demand <- use.of.REM.in.the.production.of.goods / intrinsic.demand
    
    
    
    #State (stock) variables
    ddemand.for.REM <- increase.of.the.demand.for.REM - decrease.of.demand.through.price.elasticity.of.demand - substitution.losses
    davailable.supply.of.REM <-  processing.of.REM - use.of.REM.in.the.production.of.goods
    dREM.in.goods <- use.of.REM.in.the.production.of.goods/average.lifetime.in.goods
    dextraction.capacity.under.construction <- newly.planned.extraction.capacity - commissioning.of.extraction.capacity
    dinstalled.extraction.capacity <- commissioning.of.extraction.capacity - (decommissioning.of.extraction.capacity+decommissioning.of.unprofitable.extraction.capacity)
    dintrinsic.demand <- demand.for.REM - decrease.of.demand.through.price.elasticity.of.demand - substitution.losses
    
    list(c(ddemand.for.REM,
           davailable.supply.of.REM,
           dREM.in.goods,
           dextraction.capacity.under.construction,
           dinstalled.extraction.capacity,
           dintrinsic.demand))
  })
}

parameters<-c(average.lifetime.in.goods = 15,
              real.REM.extraction= 0,
              price.elasticity.of.demand = .10,  # %
              normal.profit.margin = 0.15,
              relative.price.of.the.cheapest.substitute = 100,
              economic.growth.rate = 0.05,
              precise.construction.time.of.extraction.capacity = 8,
              average.lifetime.of.extraction.capacity = 30,
              cumulatively.extracted.REM = 4000000,
              initial.cumulatively.extracted.REM = 0) # years

InitialConditions <- c(REM.in.goods= 2000000, # ton in the year2010
                       available.supply.of.REM = 250000,
                       demand.for.REM = 100000,
                       extraction.capacity.under.construction = 60000,#tons
                       installed.extraction.capacity = 100000,
                       intrinsic.demand = 100000) 
                       
            

times <- seq(2000 , #initial time, years
             2080, #years, 
             1) #time step, years

intg.method<-c("rk4")

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

1.4. Muestra en gráficos el comportamiento del indicador ‘fraction produced of intrinsic demand’, de ‘installed extraction capacity’, de ‘scarcity price effect’, y de ‘relative price’.

plot(out, main=list("fraction.produced.of.intrinsic.demand","installed.extraction.capacity","scarcity.price.effect","relative.price"))

1.5. ¿Qué ocurre si la ‘initial extraction capacity under construction’ es 0? Compara, grafica y concluye.

library("deSolve")

REM  <- function(t, state, parameters) {
  with(as.list(c(state,parameters)), {
    #Endogenous auxiliary variables
    
    supply.demand.ratio <- available.supply.of.REM / demand.for.REM
    scarcity.price.effect <- approx(c(0,0.55,1.1,2.2,11,22),#x
                             c(100,10,1,0.75,0.5,0.2),#y
                             xout=supply.demand.ratio)$y
    
    dif.average.REM.extraction.costs <- cumulatively.extracted.REM -  initial.cumulatively.extracted.REM
    
    average.REM.extraction.costs <- approx(c(0,2000000,4000000,6000000,8000000,10000000,12000000,14000000,16000000,18000000),#x
                             c(1,2,4,8,16,32,64,128,256,512),#y
                             xout= dif.average.REM.extraction.costs)$y
    
    processing.of.REM <- real.REM.extraction
    real.REM.extraction <- ifelse(scarcity.price.effect < 1, installed.extraction.capacity*scarcity.price.effect, installed.extraction.capacity)
    
    economic.growth.rate <- if (t < 2009) {
  growth_rate <- 0.03
} else if (t == 2010) {
  growth_rate <- -0.08
} else if (t == 2011) {
  growth_rate <- 0.05
} else if (t > 2009) {
  growth_rate <- -0.10
}

    
    relative.price <- (average.REM.extraction.costs * scarcity.price.effect) *  (1 + normal.profit.margin)
    commissioning.of.extraction.capacity <- extraction.capacity.under.construction / precise.construction.time.of.extraction.capacity
      
    
    #Flow variables
    use.of.REM.in.the.production.of.goods <- max(0, demand.for.REM - available.supply.of.REM) 
    
    increase.of.the.demand.for.REM <- economic.growth.rate * demand.for.REM
    
    relative.price.of.the.previous.year <- relative.price
    
    decrease.of.demand.through.price.elasticity.of.demand <- - price.elasticity.of.demand * demand.for.REM * (((1/relative.price)-(1/relative.price.of.the.previous.year))/(1/relative.price.of.the.previous.year))
    
    substitution.losses <- max(0,ifelse(relative.price > relative.price.of.the.cheapest.substitute, demand.for.REM*(1/relative.price)*(relative.price-relative.price.of.the.cheapest.substitute), 0))
    
    desired.extraction.capacity <- demand.for.REM
    
    profitability.of.REM.extraction <- (relative.price - average.REM.extraction.costs) / average.REM.extraction.costs
    
    newly.planned.extraction.capacity <- (min(max(profitability.of.REM.extraction, 0), 1) * min(max((desired.extraction.capacity - installed.extraction.capacity), 0), installed.extraction.capacity))
    
    decommissioning.of.unprofitable.extraction.capacity <- installed.extraction.capacity * (-min(profitability.of.REM.extraction,0))
    
    decommissioning.of.extraction.capacity <- max(0,(installed.extraction.capacity / average.lifetime.of.extraction.capacity - decommissioning.of.unprofitable.extraction.capacity))
    
    fraction.produced.of.intrinsic.demand <- use.of.REM.in.the.production.of.goods / intrinsic.demand
    
    
    
    #State (stock) variables
    ddemand.for.REM <- increase.of.the.demand.for.REM - decrease.of.demand.through.price.elasticity.of.demand - substitution.losses
    davailable.supply.of.REM <-  processing.of.REM - use.of.REM.in.the.production.of.goods
    dREM.in.goods <- use.of.REM.in.the.production.of.goods/average.lifetime.in.goods
    dextraction.capacity.under.construction <- newly.planned.extraction.capacity - commissioning.of.extraction.capacity
    dinstalled.extraction.capacity <- commissioning.of.extraction.capacity - (decommissioning.of.extraction.capacity+decommissioning.of.unprofitable.extraction.capacity)
    dintrinsic.demand <- demand.for.REM - decrease.of.demand.through.price.elasticity.of.demand - substitution.losses
    
    list(c(ddemand.for.REM,
           davailable.supply.of.REM,
           dREM.in.goods,
           dextraction.capacity.under.construction,
           dinstalled.extraction.capacity,
           dintrinsic.demand))
  })
}

parameters<-c(average.lifetime.in.goods = 15,
              real.REM.extraction= 0,
              price.elasticity.of.demand = .10,  # %
              normal.profit.margin = 0.15,
              relative.price.of.the.cheapest.substitute = 100,
              economic.growth.rate = 0.05,
              precise.construction.time.of.extraction.capacity = 8,
              average.lifetime.of.extraction.capacity = 30,
              cumulatively.extracted.REM = 4000000,
              initial.cumulatively.extracted.REM = 0) # years

InitialConditions <- c(REM.in.goods= 2000000, # ton in the year2010
                       available.supply.of.REM = 250000,
                       demand.for.REM = 100000,
                       extraction.capacity.under.construction = 0,#cambiado de 6000 a 0
                       installed.extraction.capacity = 100000,
                       intrinsic.demand = 100000) 
                       
            

times <- seq(2000 , #initial time, years
             2080, #years, 
             1) #time step, years

intg.method<-c("rk4")

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

plot(out, main=list("fraction.produced.of.intrinsic.demand","installed.extraction.capacity","scarcity.price.effect","relative.price"))

Al realizar el cambio a la capacidad de extracción inicial de contstrucción a cero podemos identificar varios cambios a las otras variables claves ya que este cambio afecta directamente a la oferta disponible y a la capacidad para satisfacer la demanda intrínseca.

• Capacidad de Extracción Instalada: Cuando la capacidad inicial de extracción en construcción es cero, la capacidad instalada tarda más tiempo en incrementarse, ya que no hay capacidad inicial añadida al sistema. Esto limita la capacidad de aumentar la oferta de REM frente a una demanda creciente, afectando así otras variables del sistema.

• Efecto de la Escasez en el Precio: La falta de aumento en la capacidad de extracción instalada lleva a una mayor escasez de REM, lo que debería aumentar el efecto del precio debido a la escasez. Sin embargo, en este modelo, el efecto del precio de escasez se mantiene constante y alto, lo que indica que el precio de los REM es muy sensible a la escasez.

• Precio Relativo: En el escenario sin capacidad de extracción inicial, el precio relativo experimenta un aumento inicial antes de estabilizarse. Este aumento puede reflejar el costo inicial más alto de la extracción debido a la menor capacidad de extracción instalada y la escasez resultante.

• Fracción Producida de la Demanda Intrínseca: Se observa una disminución más marcada en la proporción de la demanda intrínseca que puede ser satisfecha en el escenario sin capacidad inicial. Esto muestra cómo la capacidad limitada de extracción afecta directamente la capacidad de satisfacer la demanda de REM.

En conclusión, reducir o eliminar la capacidad de extracción inicial en construcción tiene un impacto significativo en la capacidad del sistema para responder a la demanda de REM. Esto puede llevar a precios más altos y a una satisfacción reducida de la demanda, lo cual es crucial para las políticas de gestión de recursos naturales. Estos hallazgos pueden ayudar a los responsables de la formulación de políticas a comprender las consecuencias de las inversiones en infraestructura de extracción y la importancia de la planificación adecuada en la gestión de recursos escasos.

Problema 2: Globalización y liberalización (35 puntos)

Problema 18.19 del libro Small System Dynamics Models for Big Issues: Triple Jump towards Real-World Complexity de Erik Pruyt.

Introducción

This case is based on the Globalization model originally developed by Hartmut Bossel (2007c, Z608). In this case, 2 countries are modeled. The countries are structurally similar but differ in key values. First, country I is modeled. Adding the suffix cI to all country I variables will make copy-pasting country II much easier. That is, copy-paste (replicate with suffix I) will automatically relabel all country II variables to variables with suffix cII. Then trade barriers between the countries are suddenly lifted at globalization, at time = 10. Values in this model are relative values, e.g. all state variables are standardized to 1.

Country I

In country I, production capacity cI, initially equal to 100%, is increased through investments cI and decreased through depreciation cI. Suppose the latter variable is equal to a 5% depreciation rate cI times the production capacity cI. Model the investments cI as the product of a 10% investment rate cI, the production capacity cI, and the investment function cI, if and only if the surplus cI is negative, else as equal to 0. Assume the investment function cI is equal to 2 before globalization and to an investment factor cI afterwards. Set the investment factor cI to 2 too.

The surplus cI is defined as the supply cI minus the demand for products from cI. Assume all production capacity is used at the full 100%, i.e. supply cI is equal to the production capacity cI times a 100% production rate cI. The demand for products from cI is the sum of the domestic demand which is the percentage purchased by cI from cI times the market volume cI of 100% and the foreign demand which is the market volume cII of 100% times (1 - percentage purchased by cII from cII ). Suppose the percentage purchased by cI from cI is a function of the price ratio of domestic versus imported products in cI, i.e. price of cI products in cI / price of cII products in cI, such that the percentage purchased by cI from cI is 1 if the price ratio is 0, the percentage purchased by cI from cI is 1 if the price ratio is 0.5, the percentage purchased by cI from cI is 0.5 if the price ratio is 1, the percentage purchased by cI from cI is 0 if the price ratio is 1.5, the percentage purchased by cI from cI is 0 if the price ratio is 2, and the percentage purchased by cI from cI is 0 if the price ratio is 5.

Suppose the market in country I is so competitive that the product price cI equals the product costs times (1 + the tax rate of 20%). The product costs cI are the sum of the resource costs cI of 100% and the production costs cI which, in turn, are equal to the standards of cI times the standard factor cI of 100%. Model the standards of cI with a stock variable, initially equal to 1, which increases through progress cI and decreases through deterioration cI. Deterioration cI is proportional to the standards of cI with a deterioration rate cI of 5% per year. Progress cI could be modeled as the product of standards of cI, investment in progress cI, and (1 - product price cI /reference price cI), with a reference price cI of 5. Suppose investment in progress cI is equal to the amount of investments cI times a progress function cI. Model the progress function cI such that it is 2 until month 10 and equal to the progress factor cI thereafter. Set the progress factor cI for the moment to 2.

The price of cII products in cI is equal to the product price cII times (1 - subsidy cII + customs duties raised by cI). Assume the subsidy cI falls from 75% to 0% and the customs duties raised by cI from 50% to 0% at globalization.

Preguntas del caso, parte 1:

2.1. Para esta primera parte del caso construye el modelo de dinámica de sistemas del país I en R y simula el modelo por un periodo de 50 años.

# Carga el paquete deSolve
library(deSolve)

# Definimos los parámetros
parameters <- c(
  CUSTOMS_DUTY_c1 = 0,
  # CUSTOMS_DUTY_c2 = 0
  DEPRECIATION_RATE_c1 = 0.05,
  # DEPRECIATION_RATE_c2 = 0.05
  DETERIORATION_RATE_c1 = 0.05,
  # DETERIORATION_RATE_c2 = 0.05
  EXPORT_SUBSIDY_c1 = 0,
  # EXPORT_SUBSIDY_c2 = 0
  INVESTMENT_FACTOR_c1 = 2,
  # este se deja para modelo de solo país 1
  # INVESTMENT_FACTOR_c2 = 2
  INVESTMENT_RATE_c1 = 0.1,
  # INVESTMENT_RATE_c2 = 0.1
  MARKET_VOLUME_c1 = 1,
  MARKET_VOLUME_c2 = 1,
  # este se deja para modelo de solo país 1
  PRODUCTION_RATE_c1 = 1,
  # PRODUCTION_RATE_c2 = 1
  PROGRESS_FACTOR_c1 = 2,
  # PROGRESS_FACTOR_c2 = 2
  REFERENCE_PRICE_c1 = 5,
  # REFERENCE_PRICE_c2 = 5
  RESOURCE_COSTS_c1 = 1,
  # RESOURCE_COSTS_c2 = 1
  STANDARD_FACTOR_c1 = 1,
  # STANDARD_FACTOR_c2 = 1
  TAX_RATE_c1 = 0.2
  # TAX_RATE_c2 = 0.2
)

# Condiciones iniciales
InitialConditions <- c(
  production_capacity_c1 = 1,
  standards_c1 = 1
  # production_capacity_c2 = 0.1,
  # standards_c2 = 0.1
)

# Definimos la función del modelo
model <- function(time, state, parameters) {
  with(as.list(c(state, parameters)), {
    
    customs_c1 <- ifelse(time < 10, 0.5, CUSTOMS_DUTY_c1)
    # customs_c2 <- ifelse(time < 10, 0, CUSTOMS_DUTY_c2)
    
    production_costs_c1 <- STANDARD_FACTOR_c1 * standards_c1
    # production_costs_c2 <- STANDARD_FACTOR_c2 * standards_c2
    
    product_costs_c1 <- RESOURCE_COSTS_c1 + production_costs_c1
    # product_costs_c2 <- RESOURCE_COSTS_c2 + production_costs_c2
    
    product_price_c1 <- (1 + TAX_RATE_c1) * product_costs_c1
    # product_price_c2 <- (1 + TAX_RATE_c2) * product_costs_c2
    
    inland_price_c1 <- ifelse(time < 10, 1, product_price_c1 * (1 - EXPORT_SUBSIDY_c1 + customs_c1))
    # inland_price_c2 <- ifelse(time < 10, 1, product_price_c1 * (1 - EXPORT_SUBSIDY_c1 + customs_c2))
    
    price_ratio_domestic_vs_imported_product_c1 <- product_price_c1 / inland_price_c1
    # price_ratio_domestic_vs_imported_product_c2 <- product_price_c2 / inland_price_c2
    
    purchase_decision_c1 <- approx(c(0, 0.5, 1, 1.5, 2, 5),
                                   c(1, 1, 0.5, 0, 0, 0),
                                   xout = price_ratio_domestic_vs_imported_product_c1)$y
    
    #purchase_decision_c2 <- approx(c(0, 0.5, 1, 1.5, 2, 5),
    #                                c(1, 1, 0.5, 0, 0, 0),
    #                                xout = price_ratio_domestic_vs_imported_product_c2)$y
    
    investment_function_c1 <- ifelse(time < 10, 2, INVESTMENT_FACTOR_c1)
    # investment_function_c2 <- ifelse(time < 10, 0.2, INVESTMENT_FACTOR_c2)
    
    subsidy_c1 <- ifelse(time < 10, 0.75, EXPORT_SUBSIDY_c1)
    subsidy_c2 <- 0
    # subsidy_c2 <- ifelse(time < 10, 0, EXPORT_SUBSIDY_c2)
    
    supply_c1 <- PRODUCTION_RATE_c1 * production_capacity_c1
    # supply_c2 <- PRODUCTION_RATE_c2 * production_capacity_c2
    
    demand_c1 <- purchase_decision_c1 * MARKET_VOLUME_c1 + (1 - ifelse(time < 10, 0, 0)) * MARKET_VOLUME_c2
    # demand_c2 <- purchase_decision_c2 * MARKET_VOLUME_c2 + (1 - ifelse(time < 10, 0, purchase_decision_c1)) * MARKET_VOLUME_c1
    
    surplus_c1 <- supply_c1 - demand_c1
    # surplus_c2 <- supply_c2 - demand_c2
    
    investment_c1 <- ifelse(surplus_c1 < 0, INVESTMENT_RATE_c1 * production_capacity_c1 * investment_function_c1, 0)
    # investment_c2 <- ifelse(surplus_c2 < 0, INVESTMENT_RATE_c2 * production_capacity_c2 * investment_function_c2, 0)
    
    depreciation_c1 <- DEPRECIATION_RATE_c1 * production_capacity_c1
    # depreciation_c2 <- DEPRECIATION_RATE_c2 * production_capacity_c2
    
    deterioration_c1 <- DETERIORATION_RATE_c1 * standards_c1
    # deterioration_c2 <- DETERIORATION_RATE_c2 * standards_c2
    
    progress_function_c1 <- ifelse(time < 10, 2, PROGRESS_FACTOR_c1)
    # progress_function_c2 <- ifelse(time < 10, 1, PROGRESS_FACTOR_c2)
    
    progress_investment_c1 <- progress_function_c1 * investment_c1
    # progress_investment_c2 <- progress_function_c2 * investment_c2
    
    progress_c1 <- progress_investment_c1 * standards_c1 * (1 - product_price_c1 / REFERENCE_PRICE_c1)
    # progress_c2 <- progress_investment_c2 * standards_c2 * (1 - product_price_c2 / REFERENCE_PRICE_c2)
    
    # Variables de flujo
    
    # Variables de stock (ecuaciones diferenciales)
    d_production_capacity_c1 <- investment_c1 - depreciation_c1
    # d_production_capacity_c2 <- investment_c2 - depreciation_c2
    
    d_standards_c1 <- progress_c1 - deterioration_c1
    # d_standards_c2 <- progress_c2 - deterioration_c2
    
    # Lista de resultados
    list(c(d_production_capacity_c1, d_standards_c1))
  })
}

# Configuración temporal
times <- seq(0, 50, 0.01)

# Método de integración
intg.method <- "rk4"

# Simulación del modelo
out <- ode(
  y = InitialConditions,
  times = times,
  func = model,
  parms = parameters,
  method = intg.method
)

# Graficar los resultados
plot(out, col = c("blue"))

La capacidad de producción muestra un incremento muy rápido tras mostrarse estable los primeros 10 años, lo cual indica que la inversión y el progreso tecnológico tienen un impacto significativo. Después de este rápido aumento, la capacidad de producción se estabiliza nuevamente. Esto podría indicar que el modelo alcanza un equilibrio donde la tasa de inversión nueva se iguala con la tasa de depreciación.

Los estandares muestran un crecimiento modeerado en los primeros diez años para después crecer creca de un 60% de manera acelerada y después volver a desacelerar su crecimiento.

Country II

Assume that country II never had export subsidies, i.e. subsidy cII is zero, nor imposed customs duties. And assume that the initial standards of cII and production capacity cII were equal to 0.1. Mirroring the situation in country I, demand for products from cII is equal to the percentage purchased by cII from cII times the market volume cII plus market volume cI times (1 - purchase decision cI ). Similarly, the price of cI products in cII is equal to product price cI times (1 - subsidy cI + customs duties raised by cII). Suppose the progress function cII rises at globalization from 1 to 2 and the investment function cII from 0.2 to 2.

Preguntas del caso, parte 2:

2.2. Para esta segunda parte del caso agrega al modelo de dinámica de sistemas el país II y simula el modelo por un periodo de 50 años. La recomendación es copia y pega las variables del pais I llamadas como cI, cámbialas por CII y haz los cambios mencionados para este segundo país.

# Cargar el paquete deSolve
library(deSolve)

# Definimos los parámetros
parameters <- c(
  CUSTOMS_DUTY_c1 = 0,
  CUSTOMS_DUTY_c2 = 0,
  DEPRECIATION_RATE_c1 = 0.05,
  DEPRECIATION_RATE_c2 = 0.05,
  DETERIORATION_RATE_c1 = 0.05,
  DETERIORATION_RATE_c2 = 0.05,
  EXPORT_SUBSIDY_c1 = 0,
  EXPORT_SUBSIDY_c2 = 0,
  INVESTMENT_FACTOR_c1 = 2,
  INVESTMENT_FACTOR_c2 = 2,
  INVESTMENT_RATE_c1 = 0.1,
  INVESTMENT_RATE_c2 = 0.1,
  MARKET_VOLUME_c1 = 1,
  MARKET_VOLUME_c2 = 1,
  PRODUCTION_RATE_c1 = 1,
  PRODUCTION_RATE_c2 = 1,
  PROGRESS_FACTOR_c1 = 2,
  PROGRESS_FACTOR_c2 = 2,
  REFERENCE_PRICE_c1 = 5,
  REFERENCE_PRICE_c2 = 5,
  RESOURCE_COSTS_c1 = 1,
  RESOURCE_COSTS_c2 = 1,
  STANDARD_FACTOR_c1 = 1,
  STANDARD_FACTOR_c2 = 1,
  TAX_RATE_c1 = 0.2,
  TAX_RATE_c2 = 0.2
)

# Condiciones iniciales
InitialConditions <- c(
  production_capacity_c1 = 1,
  standards_c1 = 1,
  production_capacity_c2 = 0.1,
  standards_c2 = 0.1
)

# Definimos la función del modelo
model <- function(time, state, parameters) {
  with(as.list(c(state, parameters)), {
    customs_c1 <- ifelse(time < 10, 0.5, CUSTOMS_DUTY_c1)
    customs_c2 <- ifelse(time < 10, 0, CUSTOMS_DUTY_c2)
    
    production_costs_c1 <- STANDARD_FACTOR_c1 * standards_c1
    production_costs_c2 <- STANDARD_FACTOR_c2 * standards_c2
    
    product_costs_c1 <- RESOURCE_COSTS_c1 + production_costs_c1
    product_costs_c2 <- RESOURCE_COSTS_c2 + production_costs_c2
    
    product_price_c1 <- (1 + TAX_RATE_c1) * product_costs_c1
    product_price_c2 <- (1 + TAX_RATE_c2) * product_costs_c2
    
    inland_price_c1 <- ifelse(time<10, 1, product_price_c2 * (1 - EXPORT_SUBSIDY_c2 + customs_c1))
    inland_price_c2 <- ifelse(time<10, 1, product_price_c1 * (1 - EXPORT_SUBSIDY_c1 + customs_c2))
    
    price_ratio_domestic_vs_imported_product_c1 <- product_price_c1 / inland_price_c1
    price_ratio_domestic_vs_imported_product_c2 <- product_price_c2 / inland_price_c2
    
    purchase_decision_c1 <- approx(c(0,0.5,1,1.5,2,5), c(1,1,0.5,0,0,0), xout=price_ratio_domestic_vs_imported_product_c1)$y
    purchase_decision_c2 <- approx(c(0,0.5,1,1.5,2,5), c(1,1,0.5,0,0,0), xout=price_ratio_domestic_vs_imported_product_c2)$y
    
    investment_function_c1 <- ifelse(time < 10, 2, INVESTMENT_FACTOR_c1)
    investment_function_c2 <- ifelse(time < 10, 0.2, INVESTMENT_FACTOR_c2)
    
    subsidy_c1 <- ifelse(time < 10, 0.75, EXPORT_SUBSIDY_c1)
    subsidy_c2 <- ifelse(time < 10, 0, EXPORT_SUBSIDY_c2)
    
    supply_c1 <- PRODUCTION_RATE_c1 * production_capacity_c1
    supply_c2 <- PRODUCTION_RATE_c2 * production_capacity_c2
    
    demand_c1 <- purchase_decision_c1 * MARKET_VOLUME_c1 + (1 - ifelse(time<10,0,purchase_decision_c2)) * MARKET_VOLUME_c2
    demand_c2 <- purchase_decision_c2 * MARKET_VOLUME_c2 + (1 - ifelse(time<10,0,purchase_decision_c1)) * MARKET_VOLUME_c1
    
    surplus_c1 <- supply_c1 - demand_c1
    surplus_c2 <- supply_c2 - demand_c2
    
    investment_c1 <- ifelse(surplus_c1 < 0, INVESTMENT_RATE_c1 * production_capacity_c1 * investment_function_c1, 0)
    investment_c2 <- ifelse(surplus_c2 < 0, INVESTMENT_RATE_c2 * production_capacity_c2 * investment_function_c2, 0)
    
    depreciation_c1 <- DEPRECIATION_RATE_c1 * production_capacity_c1
    depreciation_c2 <- DEPRECIATION_RATE_c2 * production_capacity_c2
    
    deterioration_c1 <- DETERIORATION_RATE_c1 * standards_c1
    deterioration_c2 <- DETERIORATION_RATE_c2 * standards_c2
    
    progress_function_c1 <- ifelse(time < 10, 2, PROGRESS_FACTOR_c1)
    progress_function_c2 <- ifelse(time < 10, 1, PROGRESS_FACTOR_c2)
    
    progress_investment_c1 <- progress_function_c1 * investment_c1
    progress_investment_c2 <- progress_function_c2 * investment_c2
    
    progress_c1 <- progress_investment_c1 * standards_c1 * (1 - product_price_c1 / REFERENCE_PRICE_c1)
    progress_c2 <- progress_investment_c2 * standards_c2 * (1 - product_price_c2 / REFERENCE_PRICE_c2)
    
    d_production_capacity_c1 <- investment_c1 - depreciation_c1
    d_production_capacity_c2 <- investment_c2 - depreciation_c2
    
    d_standards_c1 <- progress_c1 - deterioration_c1
    d_standards_c2 <- progress_c2 - deterioration_c2
    
    list(c(d_production_capacity_c1, d_standards_c1, d_production_capacity_c2, d_standards_c2))
  })
}

# Configuración temporal
times <- seq(0, 50, 0.01)

# Método de integración
intg.method <- "rk4"

# Simulación del modelo
out <- ode(
  y = InitialConditions,
  times = times,
  func = model,
  parms = parameters,
  method = intg.method
)

# Gráfico de "production_capacity_c1" y "production_capacity_c2"
plot(out[, "time"], out[, "production_capacity_c1"], type="l", col="blue", xlab="Time", ylab="Production Capacity", ylim=range(out[, c("production_capacity_c1", "production_capacity_c2")]))
lines(out[, "time"], out[, "production_capacity_c2"], col="red")
legend("topright", legend=c("C1", "C2"), col=c("blue", "red"), lty=1)

2.3. Grafica los efectos en ‘production capacity’ y ‘standards’ para ambos países.

# Gráfico de "standards_c1" y "standards_c2"
plot(out[, "time"], out[, "standards_c1"], type="l", col="blue", xlab="Time", ylab="Standards", ylim=range(out[, c("standards_c1", "standards_c2")]))
lines(out[, "time"], out[, "standards_c2"], col="red")
legend("topright", legend=c("C1", "C2"), col=c("blue", "red"), lty=1)

2.4. En el tiempo 10, se liberaliza el comercio. ¿Qué ocurre posteriormente en ambos países? Explique.

El país 1 la capacidad de producción disminuye a partir de la liberalización del mercado dado el aumento de competitividad del país 2 que crece aceleradamente a partir de ese momento. A partir del año 25, la capacidad de producción del país 1 vuelve a incrementar, primero más aceleradamente por un corto periodo y después se estabiliza su crecimiento a una tasa menor. Por su parte, el país 2 presenta una contracción a ritmo constante a partir del años 30 aproximadamente.

En cuanto estandares, el país 1 abandona su crecimiento moderado de los primeros 10 años para comenzar a contraerse ahasta el año 25 aproximadamente, y después volver a un crecimiento acelerado hasta el año 30 y finalmente estabilizarse a una tasa de crecimiento moderada. Por su parte, el país dos comienza a tener un crecimiento constante a partir de la liberalización que a partir del año 30 comienza a estabilizarse a una misma tasa de crecimiento.

Problema 3: Colapso de las civilizaciones (35 puntos)

Problema 18.22 del libro Small System Dynamics Models for Big Issues: Triple Jump towards Real-World Complexity de Erik Pruyt.

Introducción

One of the great mysteries of human history has been the sudden collapse –around 800 AD– of one of the main centers of Mayan civilization in Central America at a time when it was apparently peaking in terms of culture, architecture and population. No one knows exactly why this society of several million people collapsed, but new research shows a gradually tightening squeeze between population and environment that may have been crucial to the fall. Tropical environments are notoriously fragile.

Research suggests that the population increased –just before the collapse– to about 200 to 500 persons per square kilometer. This population density required advanced agriculture or/and large-scale trade. Within two to four Mayan generations, the population density dropped to less than 20 persons per square kilometer (i.e. the same as 2000 years earlier). Furthermore, after the collapse, whole areas remained almost uninhabited for some thousand years. Some of the environmental changes appear to have been as long-lasting as the loss of population. Lakes that were apparently centers of settlement in the Maya era have not yet entirely recovered in terms of productivity.

Scientists of Florida State University and the University of Chicago estimate that there was an exponential growth in Mayan population during at least 1700 years in the tropical lowlands of what is now Guatemala such that the population doubled every 408 years. This trend may have caught the Maya in a strange trap. Their numbers grew at a steadily increasing pace, but, for many centuries, the growth was too slow for any single generation to see what was happening. Over the centuries, the increasing pressure on the environment may have become impossible to sustain. Yet the squeeze could have been imperceptible until the final population explosion, just before the collapse.

New estimates for the southern lowlands are based largely on a detailed survey of traces of residential structures that were built, occupied and abandoned over the centuries. The studies focuss on the region of two adjacent lakes (Lake Yaxha and Lake Sacnab) in the Peten lake district of northern Guatemala. The area was inhabited as early as 3000 years ago and the first agricultural settlements appeared there about 1000 BC. The land was largely deforested by 250 AD. Gradually-intensified agriculture and increasing settlements seem to have caused severe cumulative damage to an originally verdant environment. Essential nutrients washed away in the lakes, diminishing the fertility of agricultural land. Increases in phosphorus in the lakes from agriculture and human wastes seem to have aggravated the environmental damage.

Consideraciones

• The consumed food per person was about 400 kg per person per year.
• The Mayans deforested the area they needed to close the gap (in kg) between the food produced and the demand for food, given the fertility of the agricultural land (\(kg/km^{2}\)) and the intensity of the agriculture. Following function could be used: Deforestation <- min( gap / max(fertility of Lands, 1), forest/4) / Intensity ; Units: \(km^{2}/year\).
• Only 5% (the ‘emigration ratio’) of those that theoretically were left without food actually emigrated. The others stayed because of redistribution of the existing food.
• The stock variable fertility of agricultural land was initially 5000000 \(kg/(km^{2}*year\)). With the initial value of 5000000 \(kg/km^{2}\), 40000000 kg food could be produced with an initial 8 \(km^{2}\) of agricultural land. This is equal to the initial food demanded for a Population of 100000 persons with a consumed food per person equal to 400 kg per year. However, this initial fertility of agricultural land is reduced with a flow variable ‘fertility Losses’.
• The forest area has an initial value of 5000 \(km^{2}\) and could decrease by means of deforestation.
• Define the food Gap as demand of food - food produced.
• The intensity initially equals 1. The variable could be used to regulate the speed of deforestation and the speed of fertility losses. If intensity = 1, then all the forest needed is deforested every year.
• Use a Net natural population increase flow of the population that leads to a growth of the population from initially 100000 in 1000 BC to 2000000 inhabitants in the year 800 AD, or in other words, that leads to a doubling of the population every 408 years.
• Model the ‘fertility losses’ as: fertility of land * min(2,(agricultural land/forest)\(^{x}\))/intensity; in [\(kg/(km^{2}*year)/year\)]. This means that the fertility losses are proportional to the existing fertility of agricultural land, and to the fraction of agricultural land and Forests with an (almost) quadratic relation. Calculate the exponent \(x\) that leads to a collapse in the year 800.

Preguntas del caso:

3.1. Construye un modelo de dinámica de sistemas en R basado en el caso anterior. Calcula el valor de los exponentes y las variables de flujo omitidas de tal manera que el colapso de la población Maya ocurra en el año 800. Simula 2000 años de evolución iniciado en el año 1,000 A.C.

#install.packages("deSolve")
library("deSolve")

collapse.of.civilizations <- function(t, state, parameters) {
  with(as.list(c(state, parameters)), {
    
    # Endogenous variables
    food.produced <- fertility.of.lands * agricultural.land
    demand.for.food <- consumed.food.per.person * population
    gap <- demand.for.food - food.produced

    # Flow variables, including MIN and MAX functions to prevent unrealistic values
    deforestation <- min(gap / max(fertility.of.lands, 1), forest / 4) / intensity
    fertility.losses <- fertility.of.lands * min(2, (agricultural.land / max(forest, 1))^1.85) / intensity
    natural.population.increase <- population * natural.increase.rate
    emigration <- min((gap / consumed.food.per.person) * emigration.ratio, population * 0.05)

    # State Variables
    dagricultural.land <- deforestation
    dforest <- -deforestation
    dfertility.of.lands <- -fertility.losses
    dpopulation <- natural.population.increase - emigration

    return(list(c(dforest, dagricultural.land, dfertility.of.lands, dpopulation),
                deforestation = deforestation,
                fertility.losses = fertility.losses))
  })
}

parameters <- list(
  intensity = 1.0,
  emigration.ratio = 0.05,
  consumed.food.per.person = 400,  # kg per person per year
  natural.increase.rate = 2^(1/408) - 1
)

InitialConditions <- c(
  forest = 5000,
  agricultural.land = 8,
  fertility.of.lands = 5000000,
  population = 100000
)

times <- seq(-1000, 1000, 1)  # From 1000 BC to 2000 AD

# Run the simulation
out <- ode(y = InitialConditions, times = times, func = collapse.of.civilizations, parms = parameters, method = "rk4")

# Plotting the results
plot(out)

3.2. ¿Qué pasa con la población? Gráfica el comportamiento de esta variable, ¿Corresponde el comportamiento de tu modelo al comportamiento observado?

library(deSolve)
library(ggplot2)

# Asumiendo que 'out' contiene los resultados de la simulación ya ejecutada
results <- as.data.frame(out)
colnames(results) <- c("time", "forest", "agricultural_land", "fertility_of_lands", "population")

ggplot(results, aes(x = time, y = population)) +
  geom_line() +
  labs(title = "Población a lo largo del tiempo",
       x = "Tiempo (año)",
       y = "Población") +
  theme_minimal()

La población muestra un rápido aumento seguido de un colapso dramático, que ocurre después de que los recursos clave como el bosque se agotan y la fertilidad de la tierra disminuye significativamente. Este patrón de crecimiento seguido de un colapso repentino es consistente con lo que se espera en escenarios donde las civilizaciones enfrentan límites ambientales severos, lo cual puede ser interpretado como una simulación de lo que pudo haber contribuido al colapso de la civilización Maya alrededor del año 800 d.C.

3.3. Introduce cambios en el modelo que hagan más realista su comportamiento.

#install.packages("deSolve")
library("deSolve")

collapse.of.civilizations <- function(t, state, parameters) {
  with(as.list(c(state, parameters)), {
    
    # Endogenous variables
    food.produced <- fertility.of.lands * agricultural.land
    demand.for.food <- consumed.food.per.person * population
    gap <- demand.for.food - food.produced

    # Flow variables, including MIN and MAX functions to prevent unrealistic values
    deforestation <- min(gap / max(fertility.of.lands, 1), forest / 4) / intensity
    fertility.losses <- fertility.of.lands * min(2, (agricultural.land / max(forest, 1))^1.5) / intensity
    natural.population.increase <- population * natural.increase.rate
    emigration <- min((gap / consumed.food.per.person) * emigration.ratio, population * 0.05)

    # State Variables
    dagricultural.land <- deforestation
    dforest <- -deforestation
    dfertility.of.lands <- -fertility.losses
    dpopulation <- natural.population.increase - emigration

    return(list(c(dforest, dagricultural.land, dfertility.of.lands, dpopulation),
                deforestation = deforestation,
                fertility.losses = fertility.losses))
  })
}

parameters <- list(
  intensity = 1.2,
  emigration.ratio = 0.25,
  consumed.food.per.person = 400,  # kg per person per year
  natural.increase.rate = 2^(1/408) - 1
)

InitialConditions <- c(
  forest = 5000,
  agricultural.land = 8,
  fertility.of.lands = 5000000,
  population = 100000
)

times <- seq(-1000, 1000, 1)  # From 1000 BC to 2000 AD

# Run the simulation
out <- ode(y = InitialConditions, times = times, func = collapse.of.civilizations, parms = parameters, method = "rk4")

# Plotting the results
plot(out)

En este modelo hicimos algunos cambios, entre ellos: Intensidad de las Actividades (Intensity): Aumentamos el valor de intensity de 1.0 a 1.2. Este parámetro influye en la rapidez con la que se llevan a cabo las actividades de deforestación y las pérdidas de fertilidad, lo cual puede afectar la sostenibilidad de los recursos a largo plazo.

Ratio de Emigración (Emigration Ratio): Modificamos el ratio de emigración de 0.05 a 0.25. Este cambio representa un aumento en la proporción de la población que emigra en respuesta a la escasez de alimentos, reflejando una respuesta más drástica a las crisis de recursos.

Exponente de Pérdidas de Fertilidad (Fertility Loss Exponent): Redujimos el exponente en la fórmula de pérdidas de fertilidad de 1.85 a 1.5. Este ajuste suaviza la tasa a la que la fertilidad de las tierras disminuye en respuesta a la expansión agrícola y otros factores de estrés ambiental.

Cambios en los Resultados Observados: Bosque (Forest): La cantidad de bosque disminuye abruptamente, lo que sugiere una deforestación intensiva seguida de un cese casi total, posiblemente debido a la limitación de tierras disponibles para tal fin.

Tierra Agrícola (Agricultural Land): La tierra agrícola muestra un aumento inicial y luego se estabiliza, lo cual indica que la expansión agrícola alcanza un límite sostenible bajo las nuevas condiciones del modelo.

Fertilidad de las Tierras (Fertility of Lands): La fertilidad disminuye de manera más controlada bajo el nuevo exponente, reflejando un agotamiento gradual que se alinea mejor con los cambios en las prácticas de uso de la tierra y la deforestación.

Población (Population): La curva de población muestra un crecimiento seguido de un declive pronunciado, pero menos abrupto que en modelos anteriores, indicando un colapso poblacional que refleja mejor la interacción entre emigración y la disponibilidad decreciente de recursos. No obstante, con estas modificaciones la población en el modelo ajustado es menor en su punto máximo en comparación con el modelo anterior. Esto refleja una adaptación más realista a las condiciones limitantes ambientales y recursos decrecientes, contribuyendo a un colapso más gradual y realista de la población.

Estos cambios han contribuido a un modelo que potencialmente representa de manera más realista cómo las presiones sobre los recursos y la capacidad de carga del entorno podrían haber llevado a un colapso gradual de la civilización maya alrededor del año 800. El aumento en el ratio de emigración y la suavización en las pérdidas de fertilidad son ajustes clave que parecen haber ayudado a modelar un colapso más gradual y plausible, evitando cambios extremadamente bruscos y poco realistas en las variables del modelo.

3.4. ¿Qué tan sensible es el modelo a cambios marginales en los parámetros o a cambios estructurales? Presenta al menos dos ejemplos

El modelo muestra diferente tiene poca sensibilidad a cambios marginales por lo tanto necesita cambios significativos en los parámetros para mostrar cambios relevantes en los resultados.

Modelo Original:

#install.packages("deSolve")
library("deSolve")

collapse.of.civilizations <- function(t, state, parameters) {
  with(as.list(c(state, parameters)), {
    
    # Endogenous variables
    food.produced <- fertility.of.lands * agricultural.land
    demand.for.food <- consumed.food.per.person * population
    gap <- demand.for.food - food.produced

    # Flow variables, including MIN and MAX functions to prevent unrealistic values
    deforestation <- min(gap / max(fertility.of.lands, 1), forest / 4) / intensity
    fertility.losses <- fertility.of.lands * min(2, (agricultural.land / max(forest, 1))^1.85) / intensity
    natural.population.increase <- population * natural.increase.rate
    emigration <- min((gap / consumed.food.per.person) * emigration.ratio, population * 0.05)

    # State Variables
    dagricultural.land <- deforestation
    dforest <- -deforestation
    dfertility.of.lands <- -fertility.losses
    dpopulation <- natural.population.increase - emigration

    return(list(c(dforest, dagricultural.land, dfertility.of.lands, dpopulation),
                deforestation = deforestation,
                fertility.losses = fertility.losses))
  })
}

parameters <- list(
  intensity = 1.0,
  emigration.ratio = 0.05,
  consumed.food.per.person = 400,  # kg per person per year
  natural.increase.rate = 2^(1/408) - 1
)

InitialConditions <- c(
  forest = 5000,
  agricultural.land = 8,
  fertility.of.lands = 5000000,
  population = 100000
)

times <- seq(-1000, 1000, 1)  # From 1000 BC to 2000 AD

# Run the simulation
out <- ode(y = InitialConditions, times = times, func = collapse.of.civilizations, parms = parameters, method = "rk4")

# Plotting the results
plot(out)

Modelo con cambio de intensity = 1.2. Por si solo no muestra un gran cambio respecto al modelo original.

#install.packages("deSolve")
library("deSolve")

collapse.of.civilizations <- function(t, state, parameters) {
  with(as.list(c(state, parameters)), {
    
    # Endogenous variables
    food.produced <- fertility.of.lands * agricultural.land
    demand.for.food <- consumed.food.per.person * population
    gap <- demand.for.food - food.produced

    # Flow variables, including MIN and MAX functions to prevent unrealistic values
    deforestation <- min(gap / max(fertility.of.lands, 1), forest / 4) / intensity
    fertility.losses <- fertility.of.lands * min(2, (agricultural.land / max(forest, 1))^1.85) / intensity
    natural.population.increase <- population * natural.increase.rate
    emigration <- min((gap / consumed.food.per.person) * emigration.ratio, population * 0.05)

    # State Variables
    dagricultural.land <- deforestation
    dforest <- -deforestation
    dfertility.of.lands <- -fertility.losses
    dpopulation <- natural.population.increase - emigration

    return(list(c(dforest, dagricultural.land, dfertility.of.lands, dpopulation),
                deforestation = deforestation,
                fertility.losses = fertility.losses))
  })
}

parameters <- list(
  intensity = 1.2,
  emigration.ratio = 0.05,
  consumed.food.per.person = 400,  # kg per person per year
  natural.increase.rate = 2^(1/408) - 1
)

InitialConditions <- c(
  forest = 5000,
  agricultural.land = 8,
  fertility.of.lands = 5000000,
  population = 100000
)

times <- seq(-1000, 1000, 1)  # From 1000 BC to 2000 AD

# Run the simulation
out <- ode(y = InitialConditions, times = times, func = collapse.of.civilizations, parms = parameters, method = "rk4")

# Plotting the results
plot(out)

Modelo con cambio en el emigration.ratio por tratarse de un cambio muy drástico (de 0.05 a 0.25). En resumen, aumentar la emigration.ratio ha introducido un amortiguador en el sistema que reduce la tasa de crecimiento poblacional y modera la demanda sobre recursos naturales. Sin embargo, este cambio no ha prevenido el declive en la fertilidad del suelo, lo cual sigue siendo un problema crítico en este modelo. El comportamiento más suavizado de la población y la tierra agrícola, junto con la disminución continua en la fertilidad del suelo, reflejan un sistema que sigue enfrentando desafíos significativos para la sustentabilidad a largo plazo.

#install.packages("deSolve")
library("deSolve")

collapse.of.civilizations <- function(t, state, parameters) {
  with(as.list(c(state, parameters)), {
    
    # Endogenous variables
    food.produced <- fertility.of.lands * agricultural.land
    demand.for.food <- consumed.food.per.person * population
    gap <- demand.for.food - food.produced

    # Flow variables, including MIN and MAX functions to prevent unrealistic values
    deforestation <- min(gap / max(fertility.of.lands, 1), forest / 4) / intensity
    fertility.losses <- fertility.of.lands * min(2, (agricultural.land / max(forest, 1))^1.85) / intensity
    natural.population.increase <- population * natural.increase.rate
    emigration <- min((gap / consumed.food.per.person) * emigration.ratio, population * 0.05)

    # State Variables
    dagricultural.land <- deforestation
    dforest <- -deforestation
    dfertility.of.lands <- -fertility.losses
    dpopulation <- natural.population.increase - emigration

    return(list(c(dforest, dagricultural.land, dfertility.of.lands, dpopulation),
                deforestation = deforestation,
                fertility.losses = fertility.losses))
  })
}

parameters <- list(
  intensity = 1.0,
  emigration.ratio = 0.25,
  consumed.food.per.person = 400,  # kg per person per year
  natural.increase.rate = 2^(1/408) - 1
)

InitialConditions <- c(
  forest = 5000,
  agricultural.land = 8,
  fertility.of.lands = 5000000,
  population = 100000
)

times <- seq(-1000, 1000, 1)  # From 1000 BC to 2000 AD

# Run the simulation
out <- ode(y = InitialConditions, times = times, func = collapse.of.civilizations, parms = parameters, method = "rk4")

# Plotting the results
plot(out)

Sin embargo, si hacemos un cambio menos abrupto hacia esta variable los resultados no muestran un cambio tan notorio En este modelo cambiamos el emigration-rate de 0.05 a 0.10 y no se ve a simple vista un cambio evidente contra el modelo original.

#install.packages("deSolve")
library("deSolve")

collapse.of.civilizations <- function(t, state, parameters) {
  with(as.list(c(state, parameters)), {
    
    # Endogenous variables
    food.produced <- fertility.of.lands * agricultural.land
    demand.for.food <- consumed.food.per.person * population
    gap <- demand.for.food - food.produced

    # Flow variables, including MIN and MAX functions to prevent unrealistic values
    deforestation <- min(gap / max(fertility.of.lands, 1), forest / 4) / intensity
    fertility.losses <- fertility.of.lands * min(2, (agricultural.land / max(forest, 1))^1.85) / intensity
    natural.population.increase <- population * natural.increase.rate
    emigration <- min((gap / consumed.food.per.person) * emigration.ratio, population * 0.05)

    # State Variables
    dagricultural.land <- deforestation
    dforest <- -deforestation
    dfertility.of.lands <- -fertility.losses
    dpopulation <- natural.population.increase - emigration

    return(list(c(dforest, dagricultural.land, dfertility.of.lands, dpopulation),
                deforestation = deforestation,
                fertility.losses = fertility.losses))
  })
}

parameters <- list(
  intensity = 1.0,
  emigration.ratio = 0.10,
  consumed.food.per.person = 400,  # kg per person per year
  natural.increase.rate = 2^(1/408) - 1
)

InitialConditions <- c(
  forest = 5000,
  agricultural.land = 8,
  fertility.of.lands = 5000000,
  population = 100000
)

times <- seq(-1000, 1000, 1)  # From 1000 BC to 2000 AD

# Run the simulation
out <- ode(y = InitialConditions, times = times, func = collapse.of.civilizations, parms = parameters, method = "rk4")

# Plotting the results
plot(out)

3.5. Propón una política (i.e. preferentemente una política dinámica) que evite el colapso de la civilización Maya. Implementa esta política en el modelo y compara gráficamente el comportamiento del sistema con y sin tu política.

SISTEMA SIN POLÍTICA:

El modelo sin una intervención de política pública lleva a un colapso de la civilización alrededor del año 800 D.C. derivado del deterioro acelerado en la fertilidad de la tierra que se traduce en falta de alimentos para la población y por tanto un declive en la misma.

#install.packages("deSolve")
library("deSolve")

collapse.of.civilizations <- function(t, state, parameters) {
  with(as.list(c(state, parameters)), {
    
    # Endogenous variables
    food.produced <- fertility.of.lands * agricultural.land
    demand.for.food <- consumed.food.per.person * population
    gap <- demand.for.food - food.produced

    # Flow variables, including MIN and MAX functions to prevent unrealistic values
    deforestation <- min(gap / max(fertility.of.lands, 1), forest / 4) / intensity
    fertility.losses <- fertility.of.lands * min(2, (agricultural.land / max(forest, 1))^1.85) / intensity
    natural.population.increase <- population * natural.increase.rate
    emigration <- min((gap / consumed.food.per.person) * emigration.ratio, population * 0.05)

    # State Variables
    dagricultural.land <- deforestation
    dforest <- -deforestation
    dfertility.of.lands <- -fertility.losses
    dpopulation <- natural.population.increase - emigration

    return(list(c(dforest, dagricultural.land, dfertility.of.lands, dpopulation),
                deforestation = deforestation,
                fertility.losses = fertility.losses))
  })
}

parameters <- list(
  intensity = 1.0,
  emigration.ratio = 0.05,
  consumed.food.per.person = 400,  # kg per person per year
  natural.increase.rate = 2^(1/408) - 1
)

InitialConditions <- c(
  forest = 5000,
  agricultural.land = 8,
  fertility.of.lands = 5000000,
  population = 100000
)

times <- seq(-1000, 1000, 1)  # From 1000 BC to 2000 AD

# Run the simulation
out <- ode(y = InitialConditions, times = times, func = collapse.of.civilizations, parms = parameters, method = "rk4")

# Plotting the results
plot(out)

MODELO CON POLITICA DE MEJORA EN LA TECNOLOGIA AGRICOLA:

En este modelo exploramos la posibilidad de que a través de una política de mejora tecnológica y mejores prácticas y téncnicas agrícolas se consiga un incremento en la productividad y eficiencia en el uso de tierras agrícolas. Para lograr esto, se incorporó al modelo un factor tecnológico que incrementaba marginalmente con el paso del tiempo, así como también un ajuste reductivo a la pérdida de fertilidad de las tierras agrícolas de manera que se consiga no alterar el crecimiento poblacional y consiguiendo un mejor aprovechamiento de los recursos naturales para evitar el colapso de la civilización.

#install.packages("deSolve")
library("deSolve")

collapse.of.civilizations <- function(t, state, parameters) {
  with(as.list(c(state, parameters)), {
    
    # Ajuste de la productividad de la tierra por innovación tecnológica
    technology_factor <- 1 + 0.0005 * (t + 1000)  # Incremento gradual desde el año 0 hasta el 2000
    
    # Ajuste en la pérdida de fertilidad por mejores prácticas agrícolas
    conservation_factor <- max(0.5, 1 - 0.0001 * (t + 1000))

    # Endogenous variables
    food.produced <- fertility.of.lands * agricultural.land * technology_factor
    demand.for.food <- consumed.food.per.person * population
    gap <- demand.for.food - food.produced

    # Flow variables, including MIN and MAX functions to prevent unrealistic values
    deforestation <- min(gap / max(fertility.of.lands, 1), forest / 4) / intensity
    fertility.losses <- fertility.of.lands * min(2, (agricultural.land / max(forest, 1))^1.85) / intensity * conservation_factor
    natural.population.increase <- population * natural.increase.rate
    emigration <- min((gap / consumed.food.per.person) * emigration.ratio, population * 0.05)

    # State Variables
    dagricultural.land <- deforestation
    dforest <- -deforestation
    dfertility.of.lands <- -fertility.losses
    dpopulation <- natural.population.increase - emigration

    return(list(c(dforest, dagricultural.land, dfertility.of.lands, dpopulation),
                deforestation = deforestation,
                fertility.losses = fertility.losses))
  })
}

parameters <- list(
  intensity = 1.0,
  emigration.ratio = 0.05,
  consumed.food.per.person = 400,  # kg per person per year
  natural.increase.rate = 2^(1/408) - 1
)

InitialConditions <- c(
  forest = 5000,
  agricultural.land = 8,
  fertility.of.lands = 5000000,
  population = 100000
)

times <- seq(-1000, 1000, 1)  # From 1000 BC to 2000 AD

# Run the simulation
out <- ode(y = InitialConditions, times = times, func = collapse.of.civilizations, parms = parameters, method = "rk4")

# Plotting the results
plot(out)

3.6. Construye un diagrama de fase del modelo sin política. ¿Cuál es tu conclusión? ¿Cómo cambia este diagrama de fase con tu política?

DIAGRAMAS DE FASE SIN POLITICA: Sin utilizar ningún tipo de política, en todos los casos el crecimiento poblacional lleva al colapso cada una de las variables. Por ello, en cada una de las curvas en forma de parábola se puede observar como crece la población y los recursos naturales o la fertilidad -productividad- de estos se reduce hasta llegar prácticamente a cero.

#install.packages("deSolve")
library("deSolve")

collapse.of.civilizations <- function(t, state, parameters) {
  with(as.list(c(state, parameters)), {
    
    # Endogenous variables
    food.produced <- fertility.of.lands * agricultural.land
    demand.for.food <- consumed.food.per.person * population
    gap <- demand.for.food - food.produced
    
    # Flow variables, including MIN and MAX functions to prevent unrealistic values
    deforestation <- min(gap / max(fertility.of.lands, 1), forest / 4) / intensity
    fertility.losses <- fertility.of.lands * min(2, (agricultural.land / max(forest, 1))^1.85) / intensity
    natural.population.increase <- population * natural.increase.rate
    emigration <- min((gap / consumed.food.per.person) * emigration.ratio, population * 0.05)
    
    # State Variables
    dagricultural.land <- deforestation
    dforest <- -deforestation
    dfertility.of.lands <- -fertility.losses
    dpopulation <- natural.population.increase - emigration
    
    return(list(c(dforest, dagricultural.land, dfertility.of.lands, dpopulation),
                deforestation = deforestation,
                fertility.losses = fertility.losses))
  })
}

parameters <- list(
  intensity = 1.0,
  emigration.ratio = 0.05,
  consumed.food.per.person = 400,  # kg per person per year
  natural.increase.rate = 2^(1/408) - 1
)

InitialConditions <- c(
  forest = 5000,
  agricultural.land = 8,
  fertility.of.lands = 5000000,
  population = 100000
)

times <- seq(-1000, 1000, 1)  # From 1000 BC to 2000 AD

# Run the simulation
out <- ode(y = InitialConditions, times = times, func = collapse.of.civilizations, parms = parameters, method = "rk4")

library(deSolve)
library(ggplot2)

# Graficar el diagrama de fase agricultural_land vs. population

results <- as.data.frame(out)
colnames(results) <- c("time", "forest", "agricultural_land", "fertility_of_lands", "population")

phase_diagram <- ggplot(results, aes(x = agricultural_land, y = population)) +
  geom_path(color = "blue") +
  labs(x = "agricultural_land", y = "population") +
  theme_minimal()

print(phase_diagram)

# Graficar el diagrama de fase forest vs. population

results <- as.data.frame(out)
colnames(results) <- c("time", "forest", "agricultural_land", "fertility_of_lands", "population")

phase_diagram <- ggplot(results, aes(x = forest, y = population)) +
  geom_path(color = "blue") +
  labs(x = "forest", y = "population") +
  theme_minimal()

print(phase_diagram)

# Graficar el diagrama de fase fertility_of_land vs. population

results <- as.data.frame(out)
colnames(results) <- c("time", "forest", "agricultural_land", "fertility_of_lands", "population")

phase_diagram <- ggplot(results, aes(x = fertility_of_lands, y = population)) +
  geom_path(color = "blue") +
  labs(x = "fertility_of_lands", y = "population") +
  theme_minimal()

print(phase_diagram)

DIAGRAMAS DE FASE CON POLITICA DE MEJORA EN LA TECNOLOGIA AGRICOLA Tomando en cuenta nuestra propuesta de política de mejora en la tecnologías, prácticas y técnicas agrícolas se consigue un incremento en la población al mismo tiempo que se mejora la productividad de las tierras sin afectar los recursos naturales en su totalidad como sucede en el caso sin intervecnión. De esta manera se pudiera haber mitigado el colapso de la civilización Maya.

#install.packages("deSolve")
library("deSolve")

collapse.of.civilizations <- function(t, state, parameters) {
  with(as.list(c(state, parameters)), {
    
    # Ajuste de la productividad de la tierra por innovación tecnológica
    technology_factor <- 1 + 0.0005 * (t + 1000)  # Incremento gradual desde el año 0 hasta el 2000
    
    # Ajuste en la pérdida de fertilidad por mejores prácticas agrícolas
    conservation_factor <- max(0.5, 1 - 0.0001 * (t + 1000))

    # Endogenous variables
    food.produced <- fertility.of.lands * agricultural.land * technology_factor
    demand.for.food <- consumed.food.per.person * population
    gap <- demand.for.food - food.produced

    # Flow variables, including MIN and MAX functions to prevent unrealistic values
    deforestation <- min(gap / max(fertility.of.lands, 1), forest / 4) / intensity
    fertility.losses <- fertility.of.lands * min(2, (agricultural.land / max(forest, 1))^1.85) / intensity * conservation_factor
    natural.population.increase <- population * natural.increase.rate
    emigration <- min((gap / consumed.food.per.person) * emigration.ratio, population * 0.05)

    # State Variables
    dagricultural.land <- deforestation
    dforest <- -deforestation
    dfertility.of.lands <- -fertility.losses
    dpopulation <- natural.population.increase - emigration

    return(list(c(dforest, dagricultural.land, dfertility.of.lands, dpopulation),
                deforestation = deforestation,
                fertility.losses = fertility.losses))
  })
}

parameters <- list(
  intensity = 1.0,
  emigration.ratio = 0.05,
  consumed.food.per.person = 400,  # kg per person per year
  natural.increase.rate = 2^(1/408) - 1
)

InitialConditions <- c(
  forest = 5000,
  agricultural.land = 8,
  fertility.of.lands = 5000000,
  population = 100000
)

times <- seq(-1000, 1000, 1)  # From 1000 BC to 2000 AD

# Run the simulation
out <- ode(y = InitialConditions, times = times, func = collapse.of.civilizations, parms = parameters, method = "rk4")

library(deSolve)
library(ggplot2)

# Graficar el diagrama de fase agricultural_land vs. population

results <- as.data.frame(out)
colnames(results) <- c("time", "forest", "agricultural_land", "fertility_of_lands", "population")

phase_diagram <- ggplot(results, aes(x = agricultural_land, y = population)) +
  geom_path(color = "blue") +
  labs(x = "agricultural_land", y = "population") +
  theme_minimal()

print(phase_diagram)

# Graficar el diagrama de fase forest vs. population

results <- as.data.frame(out)
colnames(results) <- c("time", "forest", "agricultural_land", "fertility_of_lands", "population")

phase_diagram <- ggplot(results, aes(x = forest, y = population)) +
  geom_path(color = "blue") +
  labs(x = "forest", y = "population") +
  theme_minimal()

print(phase_diagram)

# Graficar el diagrama de fase fertility_of_land vs. population

results <- as.data.frame(out)
colnames(results) <- c("time", "forest", "agricultural_land", "fertility_of_lands", "population")

phase_diagram <- ggplot(results, aes(x = fertility_of_lands, y = population)) +
  geom_path(color = "blue") +
  labs(x = "fertility_of_lands", y = "population") +
  theme_minimal()

print(phase_diagram)