Imagine this

We know that another organization’s customer base (WolfCo = W) W poaches customers from M (MooseCo = M). Is it possible to model the more general question of whether and how the predation might be stalled, stymied, otherwise curtailed? This question will expand on the initial predatory interaction we saw when MooseCo meets WolfCo. In the wild world of now limited, and not so unbridled short-term profit for profit’s sake, M effectively will still do a lot of the work of customer acquisition for W while M blithely and seemingly allows the poaching campaign to continue. When and how might WolfCo just go too far?

Analogically, and in behavioral terms, are there to one team bullying another to give up its players. In a technological context are there limits to the ability of a faster, and hotter, component technology to degrade a computing platform’s performance and reliability (due to heat degradation). In polities (including any complex organization) and their statecraft (office politics gone wild) we can imagine an interaction between limits to the impact by special interests on the not-so common good.1

Given initial customer bases for each organization as well as acquisition and switching rates for each of the two organizations, we wonder

  • who leads/lags whom, by how much,

  • is the interaction reaching an equilibrium, a balance, or simple continues unabated,

  • are there any tipping points, when stock variable behavior switches direction,

  • is the process consistent with observed movements in (this case) customer acquisition and switching behavior,

  • any other interpretations?

The model we will use throughout this note is based on the Lotka (1920) - Volterra (1928) coupled interaction models more well known as predator-prey models.

Design

Our first task, as always, is to design a model, a way of systematic thinking through the interaction set-up, as a way to reaching some understanding of the interaction over time. Having reached some stable understanding will later attempt to apply a managerial (or policy) intervention.

The first figure displays notes for a design of an interaction model imagined in this note. There are six steps.

6-step framework to design a system dynamics model: limits to growth version.

  1. Pose the question for research, simulation, action, or simply discussion. Here the question continues to be how do the MooseCo and WolfCo customer bases interact given constraints on both enterprises to sustain their markets.

  2. Detail the four causes (final cause or end sought, formal cause or design, material cause or resources, and efficient cause or agents and their interventions) to answer the question. The final cause is the impact of capacity constraints on customer bases. The formal cause is a model of interactive customer base behavior over time with endogenous (within system) constraints. The material cause is the size of the customer base available to both markets as well a switching and acquisition rates. The efficient cause includes managers’ ability to change initial conditions, timing, and somehow influence switching and acquisition patterns.

  3. Develop a closed loop diagram of the interaction of stocks (e.g., customer bases here), flows ( changes in stocks ), and auxiliary components (e.g., acquisition and switch rates, initial values). See the notes for some detail. In this model we juggle three balls against the force of gravity (metaphorically!). Ball One are the Moose interacting reciprocally with Ball Two, the Wolves (of Wall Street?). Both draw from an existing Ball Three, the available size of the total market (ASM in some organizations).

  4. Develop a stock-flow-diagram to highlight the accumulate of a simulation history of flows in stocks as well as relationships among the collection of stocks and flows. Based on this diagram develop the equations for stocks, flows, and auxiliary variables. The picture on page two of the notes speaks volumes. The equations tell a similar story.

  5. Lay out a priori graphs of expected simulated behavior of stocks and flows over time. Here we depict our expectations. Due to the discounting nature of capacity constraints over time we might expect a dampening, even a smoothing, of the oscillatory behavior of the two customer bases as capacities will change as the ASM declines exponentially. But we should also watch for the acquisition and switching flows as well. They usually are quite excitable!

  6. Check that the design is a system, namely, does it possess a a) collection of relevant stocks, flows, and auxiliary variable complete enough to answer the question, (yes we have a relevant collection, but will leave doors and transoms open to add stocks, flows, auxiliary variables) b) connections to relate stocks, flows, and auxiliary variables, (several and this model is getting increasingly complicated) and c) coherent model components which support all four causes, especially the final cause, the end or good being sought (indeed, even with more complexity, we still seem to tell a reasonably interesting story of the two customer bases).

In this figure a detailed stock-flow-diagram is developed along with equations of motion of the two stocks (e.g., customer bases) through time.

Customer base interaction stock-flow-diagram.
Customer base interaction stock-flow-diagram.

We let \(t\) denote a time index with fixed intervals normalized to a size of 1, \(M(t)\) and \(W(t)\) be the stocks of MooseCo and WolfCo customers at time \(t\) with initial (\(t=0\)) values \(M(0)\) and \(W(0)\), and with subscripts \(m\) and \(w\) for MooseCo and Wolfco we denote customer acquisition rates \(a_m\) and \(a_w\) and switching rates \(s_m\) and \(s_w\), all between 0 and 1.

\[ M(t) = M(t-1) + a_m (1 - M(t-1)/K_M)M(t-1) - a_m M(t-1)W(t-1), \] and

\[ W(t) = W(t-1) + s_m a_w (1-W(t-1)/K_W)M(t-1)W(t-1) - s_w W(t-1). \] These equations determine the simultaneous interactions through the coupling moderation terms in \(M(t-1)W(t-1)\).

When we add a total market available base of customers \(C(t)\) we will add three more equations. One for a time dependent allocation rule, \(\phi(t)\), based on the ratio of MooseCo to total market (MooseCo plus WolfCo) customers.

\[ C(t) = C(t-1) + [a_m (1 - M(t-1)/K_M)M(t-1) s_m a_w (1-W(t-1)/K_W)M(t-1)W(t-1)]C(t-1) \] Complicated! The allocation rule is a little easier on the eyes.

\[ \phi(t) = \frac{M(t)}{M(t) + W(t)} \]

To determine the MooseCo time \(t\) capacity \(K_M(t)\) we incorporate the customer stock at simulation time \(t\), \(C(t)\). Similarly, we calculate the WolfCo time \(t\) capacity \(K_W(t)\).

\[ K_M(t) = \phi(t) C(t) \] \[ K_W(t) = (1-\phi(t)) C(t) \]

We can well ask how many exogenous auxiliary variables remain in the system. And then are these variables amenable to managerial influence.

Simulate

We continue to feature Vensim and R simulations. The model we develop here is still based on McElreath (2020), Chapter 16, p. 541-550, much modified by the capacity discounts to the acquisition of \(M\) or \(W\) customers.

Vensim model

Stocks and flows

This initial model renders the influences exogenous capacity constraints. These constraints are by construction here also not interactive. Stock-flow-model In the next model, we endogenize capacity constraints in a dynamic way by relating allocation of available customers through the time \(t\) levels of \(M(t)\) and \(W(t)\).

Endogenized Capacity Limits to Growth.
Endogenized Capacity Limits to Growth.

A very different set of behaviors befalls our model.

Equations

(01)    aCustomersInitial=
        100
    Units: **undefined** [35,150]
    
(02)    aMooseAcquisitionRate=
        0.5
    Units: **undefined**
    
(03)    aMooseAllocationFraction=
        sMoose/(sMoose+sWolf)
    Units: **undefined**
    
(04)    aMooseAvailability=
        1-sMoose/aMooseCapacity
    Units: **undefined**
    
(05)    aMooseCapacity=
        aMooseAllocationFraction * sCustomers
    Units: **undefined**
    
(06)    aMooseInitial=
        30
    Units: **undefined**
    
(07)    aMooseSwitchingRate=
        0.05
    Units: **undefined**
    
(08)    aWolfAcquisitionRate=
        0.5
    Units: **undefined**
    
(09)    aWolfAvailability=
        1-sWolf/aWolfCapacity
    Units: **undefined**
    
(10)    aWolfCapacity=
        (1-aMooseAllocationFraction)*sCustomers
    Units: **undefined**
    
(11)    aWolfInitial=
        4
    Units: **undefined**
    
(12)    aWolfSwitchingRate=
        0.5
    Units: **undefined**
    
(13)    fCustomer=
        aMooseAcquisitionRate+aWolfAcquisitionRate
    Units: **undefined**
    
(14)    FINAL TIME  = 24
    Units: Month
    The final time for the simulation.

(15)    fMooseAcquisition=
        aMooseAcquisitionRate * aMooseAvailability * sMoose
    Units: **undefined**
    
(16)    fMooseSwitching=
        aMooseSwitchingRate * sMoose * sWolf
    Units: **undefined**
    
(17)    fWolfAcquisition=
        aMooseSwitchingRate * aWolfAcquisitionRate * aWolfAvailability * sMoose *
     sWolf
    Units: **undefined**
    
(18)    fWolfSwitching=
        aWolfSwitchingRate * sWolf
    Units: **undefined**
    
(19)    INITIAL TIME  = 0
    Units: Month
    The initial time for the simulation.

(20)    SAVEPER  = 
            TIME STEP
    Units: Month [0,?]
    The frequency with which output is stored.

(21)    sCustomers= INTEG (
        -fCustomer,
            aCustomersInitial)
    Units: **undefined**
    
(22)    sMoose= INTEG (
        fMooseAcquisition-fMooseSwitching,
            aMooseInitial)
    Units: **undefined**
    
(23)    sWolf= INTEG (
        fWolfAcquisition-fWolfSwitching,
            aWolfInitial)
    Units: **undefined**
    
(24)    TIME STEP  = 0.0625
    Units: Month [0,?]
    The time step for the simulation.

Both time series (top) and phase diagram (bottom) show repeating cyclic (oscillatory) behavior, one leading / lagging the other. But now with available customers the large amplitudes and high frequencies are smoothed. It is time to work these models for further understanding.

R model

Build simulation function

Using the equations we construct a simulation function in R.

#
library( tidyverse )
#
sim_Wolf_Moose_limits <- function( n_steps, initial, gamma, dt = 0.002 ){
  W <- rep( NA, n_steps )
  M <- rep( NA, n_steps )
  W[1] <- initial[1]
  M[1] <- initial[2]
  for ( t in 2:n_steps ){
    M[t] <- M[t-1] + (gamma[1]*(1-M[t-1]/gamma[5])*M[t-1] - gamma[2]*M[t-1]*W[t-1]) * dt
    W[t] <- W[t-1] + (gamma[3]*gamma[2] * (1-W[t-1]/gamma[6]) * M[t-1] * W[t-1] - gamma[4]*W[t-1]) * dt
  }
  return( tibble(
    W = W,
    M = M
  ))
}

We introduced a time step interval size \(dt\) which will simulate \(n-steps\) values of \(M(t)\) and \(W(t)\).

Simulation data

Having defined a function, we proceed to make some simulation data for the function including a vector of acquisition and switching rates, and a vector of initial customer base values.

gamma <- c( 0.50, 0.05, 0.50, 0.50, 150, 30 )
# c( moose acquisition rate, moose switching, wolf acquistion, wolf switching, moose capacity, wolf capacity)
initial_W_M <- c( 4, 30 ) # 4 WolfCo and 30 MooseCo customers

Run simulations (until morale improves!)

In this model we will simulate a bit more carefully. We run the model with time step dt=1/30, a one day fraction of a 30 day month, and n_steps=24/(1/30)=24*30=720 days across 24 months.

days_in_month <- 30
dt_step <- 1 / days_in_month
n_steps <- 24 * days_in_month # just to put a point on it
z <- sim_Wolf_Moose_limits( n_steps, initial_W_M, gamma, dt_step )
str( z )
## tibble [720 × 2] (S3: tbl_df/tbl/data.frame)
##  $ W: num [1:720] 4 4.02 4.04 4.06 4.08 ...
##  $ M: num [1:720] 30 30.2 30.4 30.6 30.8 ...

We should plot our handiwork next.

Plot it!

day <- c(1:n_steps, 1:n_steps)
customers <- c(z$W, z$M)
group <- c(rep("WolfCo", n_steps), rep("MooseCo", n_steps))
data_plot <- tibble(
  day=day,
  customers=customers,
  group=group
)
library(plotly)
p <- data_plot |>
  ggplot( aes( x=day, y=customers, color = group) ) +
  geom_line() +
  facet_wrap( ~ group, nrow=2, scales="free" ) +
  xlab( "days" ) + ylab( "customers" ) +
  ggtitle( "Curtailing MooseCo and WolfCo customers")
ggplotly( p )

The facet_wrap layer produces over and under plots of the MooseCo and WolfCo simulated customer bases. They have very similar patterns but peaks as we would expect: first MooseCo, then WolfCo. What causes the oscillations? What causes their dampening?

Integrating Vensim and R

In this step we replicate the model, but this time read the Vensim model into R using the readsdr package. We will have to install the development package from the Github site with another package called devtools. With all of that in hand we then attempt to solve the model using the deSolve R package.

library(readsdr)
library(tidyverse)
library(plotly)
# This path was copied from the file directory properties
#  C:\Users\wgfoo\OneDrive\Documents\00000-strategic-decisions\0-models\site\models
# `\` is a reserved token called an escape. In order to read only a `\` we have to escape the escape token.
# Thus `\\`. We can automate this with grep(); 
# for another time and place
# 
model_path <- "C:\\Users\\wgfoo\\OneDrive\\Documents\\00000-strategic-decisions\\0-models\\site\\models\\wolf-moose-limits.xmile"
mdl <- read_xmile( model_path )
#
stocks_df <- sd_stocks(mdl)
stocks_df
##         name init_value
## 1 sCustomers        150
## 2     sMoose         30
## 3      sWolf          4
#
# mutate() will effectively put the data into decimal format
constants_df <- sd_constants(mdl) |> 
  mutate(value = format(value, scientific = FALSE))
constants_df
##                    name  value
## 1     aCustomersInitial 150.00
## 2 aMooseAcquisitionRate   0.50
## 3         aMooseInitial  30.00
## 4   aMooseSwitchingRate   0.05
## 5  aWolfAcquisitionRate   0.50
## 6          aWolfInitial   4.00
## 7    aWolfSwitchingRate   0.50
#
equations_df <- mdl$deSolve_components$func
equations_df
## function (time, stocks, auxs) 
## with(as.list(c(stocks, auxs)), {
##     aMooseAllocationFraction <- sMoose/(sMoose + sWolf)
##     aMooseCapacity <- aMooseAllocationFraction * sCustomers
##     aWolfCapacity <- (1 - aMooseAllocationFraction) * sCustomers
##     fMooseSwitching <- aMooseSwitchingRate * sMoose * sWolf
##     fWolfSwitching <- aWolfSwitchingRate * sWolf
##     aMooseAvailability <- 1 - sMoose/aMooseCapacity
##     aWolfAvailability <- 1 - sWolf/aWolfCapacity
##     fMooseAcquisition <- aMooseAcquisitionRate * aMooseAvailability * 
##         sMoose
##     fWolfAcquisition <- aMooseSwitchingRate * aWolfAcquisitionRate * 
##         aWolfAvailability * sMoose * sWolf
##     fCustomer <- fMooseAcquisition + fWolfAcquisition
##     d_sCustomers_dt <- -fCustomer
##     d_sMoose_dt <- fMooseAcquisition - fMooseSwitching
##     d_sWolf_dt <- fWolfAcquisition - fWolfSwitching
##     return(list(c(d_sCustomers_dt, d_sMoose_dt, d_sWolf_dt), 
##         aMooseAllocationFraction = aMooseAllocationFraction, 
##         aMooseCapacity = aMooseCapacity, aWolfCapacity = aWolfCapacity, 
##         fMooseSwitching = fMooseSwitching, fWolfSwitching = fWolfSwitching, 
##         aMooseAvailability = aMooseAvailability, aWolfAvailability = aWolfAvailability, 
##         fMooseAcquisition = fMooseAcquisition, fWolfAcquisition = fWolfAcquisition, 
##         fCustomer = fCustomer, aCustomersInitial = aCustomersInitial, 
##         aMooseAcquisitionRate = aMooseAcquisitionRate, aMooseInitial = aMooseInitial, 
##         aMooseSwitchingRate = aMooseSwitchingRate, aWolfAcquisitionRate = aWolfAcquisitionRate, 
##         aWolfInitial = aWolfInitial, aWolfSwitchingRate = aWolfSwitchingRate))
## })
## <environment: 0x0000022f3d224ff0>

That last read out is the model in all its glory. We should compare the equations with the Vensim equations. Where are the Vensim INTEG() calls? We an find them in the d_X_dt equations.

A simpler way to simulate?

We will again read the model from XMILE and then transform the contents into a form readable by deSolve using the xmile_to_deSolve() function. And as the world continues to revolve around the sun, the development version of readsdr has the sd_simulate() function to pack nearly all of what we just accomplished into a one-stop-does-all workflow. We can test drive this marvel right here.

# This path was copied from the file directory properties
# `\` is a reserved token called an escape. In order to read only a `\` we have to escape the escape token.
# Thus `\\`. We can automate this with grep(); 
# for another time and place
# 
model_path <- "C:\\Users\\wgfoo\\OneDrive\\Documents\\00000-strategic-decisions\\0-models\\site\\models\\wolf-moose-limits.xmile"
mdl <- read_xmile( model_path )
deSolve_components <- mdl$deSolve_components
dt <- 1/30
simtime <- seq(deSolve_components$sim_params$start,
               deSolve_components$sim_params$stop,
               dt)
# simulate model
result_df_endo <- ode(y = deSolve_components$stocks,
                  times  = simtime,
                  func   = deSolve_components$func,
                  parms  = deSolve_components$consts, 
                  method = "rk4")
#ds_inputs <- xmile_to_deSolve( model_path )
#ds_inputs$consts[["aWolfInitial"]] <- 4
#result_df_2 <- sd_simulate(ds_inputs, 
#            0, 
#            24, 
#            0.03333333, 
#            "rk4")
#summary( result_df_2 )

Now to use the tools we have already developed (maybe a function is in our future!) we plot our results.

z <- as.data.frame(result_df_endo)
n_steps <- length( z$time )
day <- c(1:n_steps, 1:n_steps, 1:n_steps, 1:n_steps, 1:n_steps, 1:n_steps)
customers <- c(z$sCustomers, z$fCustomer, z$sMoose, z$fMooseAcquisition, z$sWolf, z$fWolfAcquisition)
group <- c( rep("Market Base", n_steps), rep("Market Flow", n_steps), rep("MooseCo Base", n_steps), rep("MooseCo Flow", n_steps), rep("WolfCo Base", n_steps), rep("WolfCo Flow", n_steps))
data_plot_ODE <- tibble(
  day=day,
  customers=customers,
  group=group
)
library(plotly)
p <- data_plot_ODE |>
  ggplot( aes( x=day, y=customers, color=group ) ) +
  geom_line( ) +
  facet_wrap( ~group, nrow=3, scales="free") +
  xlab( "days" ) + ylab( "customers" ) +
  ggtitle( "Curtailing MooseCo and WolfCo customers")
ggplotly( p )

We might contemplate a further look into incidences of customer base per day. To do this we need only perform a little bit of arithmetic and a slight change in group names to get the idea of a net flow across our desks.

n_steps <- length( z$time )
day <- c(1:n_steps, 1:n_steps, 1:n_steps, 1:n_steps, 1:n_steps, 1:n_steps)
customers <- c(z$sCustomers, z$fCustomer, 
               z$sMoose, z$fMooseAcquisition - z$fMooseSwitching, z$sWolf, z$fWolfAcquisition - z$fWolfSwitching )
# Yes, the net flow is the difference between an
# acquisition flow and a switching flow
group <- c( rep("Market Base", n_steps), rep("Market Flow", n_steps), rep("MooseCo Base", n_steps), rep("MooseCo Net Flow", n_steps), rep("WolfCo Base", n_steps), rep("WolfCo Net Flow", n_steps))
data_plot_ODE <- tibble(
  day=day,
  customers=customers,
  group=group
)
#library(plotly)
# No changes here except perhaps for the title?
p <- data_plot_ODE |>
  ggplot( aes( x=day, y=customers, color=group ) ) +
  geom_line( ) +
  facet_wrap( ~group, nrow=3, scales="free") +
  xlab( "days" ) + ylab( "customers" ) +
  ggtitle( "MooseCo and WolfCo Incidents")
ggplotly( p )

The net flows for the two customer bases seem deeper and even negative at times. What would negative mean? In this context there is more switching out of a customer base and acquiring. The net flow seems to reach a near zero level for both customer bases. What sense does this make!? There is plenty more market available. Has time run out for whatever MooseCo and WolfCo sell?

Integration under the hood

The rk4 integration option we deployed with the sd_simulate() function is the 4th Order Runge-Kutta approximation technique introduced by Runge (1895) to integrate a first order differential equation. We can suppose this simple scalar first order homogeneous differential equation as an example.

\[ x'(t) = ax(t) \] with initial, time \(t=0\), condition

\[ x(0) = x_0. \]

For WolfCo customer acquisition \(a=0.50\) and \(x(0)=4\).

a <- 0.5                # WolfCo acquisition rate
x0 <-  4                # Initial Condition
h <- 1/30               # Time step, one day at a time
t <-  seq( 0, 24, h)    # t goes from 0 to 24 months 
xanalytical = 4*exp(0.50*t)    # Analytical solution 
# In general we won't know this, but it helps to calibrate our expectations! And this is not exact,
# since exp(x) is an approximation as well.
xstar <-  rep( 0, length(t) ); # initiate array (good coding practice to reserve some precious memory)
xstar[1] <-  x0         # Initial condition
# starts off our series of solutions at t=0
#
for (i in 1:(length(t)-1)) {
  k1 <-  a*xstar[i]  
  # First of 3 RK slopes approximates derivative
  x1 <-  xstar[i]+k1*h/2      
  # Intermediate value (using k1)
  k2 <-  a*x1        
  # A slightly better approximate derivative
  x2 <-  xstar[i]+k2*h/2      
  # Honing down an even better approximation
  k3 <-  a*x2         
  x3 <-  xstar[i]+k3*h
  # End point value 
  k4 = a*x3
  # Approximate derivative again, at end point value
  xstar[i+1] = xstar[i] + (k1+2*k2+2*k3+k4)*h/6 
  # the 4th (order) solution
  # Looks a little like a Pascal sequence: 1, 2, 2, 1
}
month <- c(t, t, t)
difference <- xanalytical - xstar
x <- c( xstar, xanalytical, difference )
method <- c( rep( "RK4", length(t)), rep( "analytical", length(t)), rep( "difference", length(t)))
data_plot_rk4 <- data.frame( month=month,
                             x=x,
                             method=method)
p <- data_plot_rk4 |>
  ggplot( aes( x=month, y=x, color=method) ) +
  geom_line() +
  facet_wrap( ~method, nrow=3, scales = "free" ) +
  xlab( "month t") + ylab( "x(t)") +
  ggtitle( "Approximate RK4 solution of x'(t) = 0.50 x(t), x(0) = 4.")
p

The so-called exact formula is simply and mathematically called analytic, that is, defined on a formal power series. There is a difference between the \(exp()\) power series expression and the rk4, effectively a difference in approximations! There is then no basis in the formulations to choose one or the other.

Generative Participative Transformation (real GPT))

This is not a large language model! But we are generating data through the interaction (participation?) of two entities which seems to transform the market.

There are two components to the generative modeling process. One is the “science” contained in the system of ordinary differential equations we just loaded into R from Vensim and then simulated with visualization. This is the process through which we formulate hypotheses about the data. The key factors include the constant parameters and the structure of the model itself.

So far we have peered into three structures, four if we count no interaction. They include simple interaction, constant capacity interaction, and endogenous capacity interaction with an available market. The hypotheses of surrounding model structure can, as we have seen, be quite complex.

The second component is the observational model, which some refer to as the measurement model. We hypothesize that this model will accurately track observed data. For us the data will be the WolfCo acquisition flow. Statistically is just an integer count of customers WolfCo acquires at each day in the simulation history. Integer counts can be modeled by the Poisson family of distributions well enough to illustrate our approach. We will eventually use the Negative Binomial distribution to account for a massive inability of the Poisson distribution to account for the variance of counts in reasonable settings like predation or customer bases. More on this timely and very sensitive topic in future episodes.

For now we will generate a set of life-like trajectories of customer base data by simply sampling from two of the 7 constants, now called parameters and parm (not eggplant parm!) in the script below. The vibrations of these two prior parameters will work their way through the labyrinth of our interactive system of ordinary differential equations to produce posterior outputs. We do not know what distribution these outputs produce, but likely they are some sort of variation on the Poisson theme.

In this exercise we will replicate as best we can the results of the Vensim model.

library(deSolve)
#
model_path <- "wolf-moose-limits.xmile"
mdl <- read_xmile( model_path )

# Pull out the ODE model from mdl
deSolve_components <- mdl$deSolve_components
# Set up initial ODE standing parameters
# widen dt to every 6 day week in a 30 day month
# in 24 months this will yield 120 6-day periods
dt <- 6/30
# seq builds a vector of times as 6 day fractions of a 
# 30 day month (dt) from 0 (start) to month 24 (stop)
simtime <- seq(deSolve_components$sim_params$start,
               deSolve_components$sim_params$stop,
               dt)
              #deSolve_components$sim_params$dt)
# parameters we will "vibrate"
# parm_names <- list("aWolfAcquisitionRate", "aMooseAcquisitionRate") for future parsing
n_steps <- deSolve_components$sim_params$stop / dt 
# 120
# initialize generative variables
sWolf_gen <- 0*(1:n_steps)
sMoose_gen <- sWolf_gen
sCustomers_gen <- sWolf_gen
# draw samples, run ODE model, save results n_steps times
for (i in 1:n_steps) {
# sample from uniform distribution == very diffuse!
# nearly equivalent to a sensitivity analysis, but
# with a probabilistic structure now
deSolve_components$consts[["aWolfAcquisitionRate"]] <- runif(1, min=0.25, max=0.75)
deSolve_components$consts[["aMooseAcquisitionRate"]] <- runif(1, min=0.25, max=0.75)
#deSolve_components$consts[["aWolfSwitchingRate"]] <- runif(1, min=0.25, max=0.75)
#deSolve_components$consts[["aMooseSwitchingRate"]] <- runif(1, min=0.025, max=0.075)
# solve
out_deSrk4 <- ode(y = deSolve_components$stocks,
                  times  = simtime,
                  func   = deSolve_components$func,
                  parms  = deSolve_components$consts, 
                  method = "rk4")
sWolf_gen[i] <- out_deSrk4[i, "sWolf"]
sMoose_gen[i] <- out_deSrk4[i, "sMoose"]
sCustomers_gen[i] <- out_deSrk4[i, "sCustomers"]
} # Thus endeth the generation
# wrangle plot data
results_gen <- data.frame(time = 1:n_steps,
                          sWolf = sWolf_gen,
                          sMoose = sMoose_gen,
                          sCustomers = sCustomers_gen
                          )
summary( results_gen )
##       time            sWolf              sMoose        sCustomers    
##  Min.   :  1.00   Min.   :0.000225   Min.   :13.94   Min.   : 18.95  
##  1st Qu.: 30.75   1st Qu.:0.066289   1st Qu.:25.74   1st Qu.: 43.06  
##  Median : 60.50   Median :1.043298   Median :37.79   Median : 57.76  
##  Mean   : 60.50   Mean   :2.143535   Mean   :38.35   Mean   : 63.59  
##  3rd Qu.: 90.25   3rd Qu.:3.909638   3rd Qu.:50.77   3rd Qu.: 74.93  
##  Max.   :120.00   Max.   :9.351079   Max.   :66.10   Max.   :150.00
# write.csv(results_gen, "results_gen.csv")

We might do well to recall that the times are six day fractions of a 30 day month, 24 months over. Now to plot our handiwork. we will use a long format. This stacks columns of data the sorts each component in the stack with an identifying key, here called group. We note our penchant for using some good coding practices here. For example we will not hard code the number of time steps n_steps, instead letting the data tell us how many there are. If we do the hard work of setting up our plotting data in this manner, the ggplot2 functions will in turn do the heavy lifting of what to put where and with what visual effects, including scaling.

z <- results_gen
n_steps <- length( z$time )
day <- c(1:n_steps, 1:n_steps)
customers <- c(z$sWolf, z$sMoose)
group <- c( rep("WolfCo Base", n_steps), rep("MooseCo Base", n_steps))
data_plot_gen <- tibble(
  day=day,
  customers=customers,
  group=group
)
library(plotly)
p <- data_plot_gen |>
  ggplot( aes( x=day, y=customers, color=group ) ) +
  geom_line( ) +
  facet_wrap( ~group, nrow=2, scales="free") +
  xlab( "6 day period" ) + ylab( "customers" ) +
  ggtitle( "Curtailing MooseCo and WolfCo customers")
ggplotly( p )

This is one example of what the actual customer data might look like. We used very uninformative priors to sample acquisition rates. There are other constants which might be sampled as well. We should look at what amounts to the (nearly, as yet to be normalized) posterior distributions of WolfCo and MooseCo customer bases. We can accomplish part of this task with a 2 dimensional density diagram to display the contours of the distribution of the scatter plot interaction data.

Instead of the long format grouping of data, we will rely on the initial short format data frame (columns of separate data) in results_gen. Here goes.

p <- results_gen |>
  ggplot( aes( x=sWolf, y=sMoose)) +
  geom_point( color = "blue" ) +
  geom_quantile(quantiles = c(0.10, 0.90)) + 
  geom_quantile(quantiles = 0.5, 
                linetype = "longdash") +
  geom_density_2d(colour = "red") +
  xlab("WolfCo Customers") + 
  ylab("MooseCo Customers") +
  ggtitle("WolfCo-MooseCo Customer Base Interaction") 
ggplotly( p )
## Smoothing formula not specified. Using: y ~ x
## Smoothing formula not specified. Using: y ~ x

Very dramatic clusters appear in this relationship at low levels of WolfCo predation of MooseCo customers. There is much to be learned from this diagram. We can begin with the credibility intervals set at 10% and 90% (cumulative probability with a purely probabilistic interpretation. Here we can say that consistent with the data of this generative model there is a 80% probability of a negative interactive relationship. The lines represent the quantile regressions of MooseCo versus WolfCo. The median 50% is the median quantile regression with 50% probability that generated observations are above (below) the line.

What have we accomplished?

Quite a bit!

  • Upgraded the Vensim model for capacity limits to growth;

  • Further upgraded the model to endogenize the capacities by introducing a market from which MooseCo and WolfCo can draw customers;

  • We then solved the model directly in R, repurposing our previous work and producing an interactive visualization using ggplot2 and plotly; we happened to use an integration method called euler;

  • We learned how to read the Vensim model into R using the readsdr package and peer into its contents;

  • From there we managed to solve the model using the deSolve package and the rk4, 4th order Runge-Kutta approximate integration technique, which we also investigated a little bit;

  • Last and the penultimate goal of all of these exercises we generated data by rerunning the model each time taking a sample of prior unobserved data from the acquisition rates;

  • Let’s not forget the 2 dimension scatter plot of generated data using the geom_density_2d() and geom_quantile() features from ggplot2.

An interactive flexdashboard app is available to play with various ranges of parameter and initial values: https://wgfoote.shinyapps.io/moose-meets-wolf-live/#section-endogenize-capacity

To do

These are items to do in no particular order.

  • Unpack all we just went through, replicate the code, play by changing parameter values and other assumptions of the model.

  • Take 3 deep breaths!

  • Review and start to execute the requirements of the assignment.

  • Review your basic statistics notes and revisit the Poisson distribution, joint probability, formation of hypotheses, product rule for both-and probability statements, sum rule for either-or probability statements.

  • When you get a chance look up the negative binomial, the gamma, the beta, the exponential, and the asymmetric Laplace distributions for reference.

The Assignment

For the Public Benefit Corporation example from Week 1,

  1. Modify the Vensim model for Week 2 to comply with your example. Describe the modifications.

  2. Read the Vensim model into R and generate data from the model adding sampling from aMooseSwitchingRate and aWolfSwitchingRate. Visualize the generated results.

  3. Interpret the generated joint distribution of the interaction of the WolfCo and MooseCo customer bases.

We will save polynumbers ( a generalization for polynomials ) for another time and space.

\[ \pi=\left\vert \overline{% \begin{array} [c]{c}% u_{0}\\ u_{1}\\ \vdots\\ u_{k}% \end{array} }\right. \]

References

Lotka, Alfred J. 1920. “Analytical Note on Certain Rhythmic Relations in Organic Systems.” Proceedings of the National Academy of Sciences 6 (7): 410–15.
McElreath, Richard. 2020. Statistical Rethinking: A Bayesian Course with Examples in R and Stan. Second edition. CRC Press. https://xcelab.net/rm/statistical-rethinking/.
Runge, Carl. 1895. “Über Die Numerische Auflösung von Differentialgleichungen.” Mathematische Annalen 46 (2): 167–78.
Volterra, Vito. 1928. “Variations and Fluctuations of the Number of Individuals in Animal Species Living Together.” ICES Journal of Marine Science 3 (1): 3–51.

  1. A bit more blunt is this quote from James K. Galbraith: “The predator-prey model explains some things that other models cannot: in particular, cycles of prosperity and depression. Growth among the prey stimulates predation. The two populations grow together at first, but when the balance of power shifts toward the predators (through rising interest rates, utility rates, oil prices, or embezzlement), both can crash abruptly. When they do, it takes a long time for either to recover.” James K. Galbraith, (retrieved 2024-03-08 from https://www.motherjones.com/politics/2006/05/predator-state/) Perhaps a limits to growth model might shave some of the extremes of boom and bust worldviews.↩︎