Making Connections

There are times when it is helpful to calibrate, or fit, your model to historical data. This capability is not built into the iThink/STELLA program, but it is possible to interface to external programs to accomplish this task. One generally available program to calibrate models is PEST, available freely from www.pesthomepage.org. In this blog post, I will demonstrate how to calibrate a simple STELLA model using PEST on Windows. Note that this method relies on the Windows command line interface added in version 9.1.2 and will not work on the Macintosh. The export to comma-separated value (CSV) file feature, added in version 9.1.2, is also used.

The model and all files associated with its calibration are available by clicking here.

The Model

The model being used is the simple SIR model first presented in my blog post Limits to Growth. The model is shown again below. There are two parameters: infection rate and recovery rate. Technically, the initial value for the Susceptible stock is also a parameter. However, since this is a conserved system, we can make an excellent guess as to its value and do not need to calibrate it.

The Data Set

We will calibrate this model to two data sets. The first is the number of weekly deaths caused by the Hong Kong flu in New York City over the winter of 1968-1969 (below).

The second is the number of weekly deaths per thousand people in the UK due to the Spanish flu (H1N1) in the winter of 1918-1919 (shown later).

In both cases, I am using the number of deaths as a proxy for the number of people infected, which we do not know. This is reasonable because the number of deaths is directly proportional to the number of infected individuals. If we knew the constant of proportionality, we could multiply the deaths by this constant to get the number of people infected.

Read more…

The Shifting the Burden Systems Archetype shows how attacking symptoms, rather than identifying and fixing fundamental problems, can lead to a further dependence on symptomatic solutions.  This Systems Archetype was formally identified in Appendix 2 of The Fifth Discipline by Peter Senge (1990).  The Causal Loop Diagram (CLD) is shown below.

When a problem symptom appears, two options present themselves:  1) apply a short-term fix to the symptom, or 2) identify and apply a longer-term fix to the fundamental issue.  The second option is less attractive because it involves a greater time delay and probably additional cost before the problem symptom is relieved.  However, applying a short-term fix, as a result of relieving the problem symptoms sooner, reduces the desire to identify and apply a more permanent fix.  Often the short-term fix also induces a secondary unintended side-effect that further undermines any efforts to apply a long-term fix.  Note that the short-term fix only relieves the symptoms, it does not fix the problem.  Thus, the symptoms will eventually re-appear and have to be addressed again.

Classic examples of shifting the burden include:

• Making up lost time for homework by not sleeping (and then controlling lack of sleep with stimulants)
• Borrowing money to cover uncontrolled spending
• Feeling better through the use of drugs (dependency is the unintended side-effect)
• Taking pain relievers to address chronic pain rather than visiting your doctor to try to address the underlying problem
• Improving current sales by focusing on selling more product to existing customers rather than expanding the customer base
• Improving current sales by cannibalizing future sales through deep discounts
• Firefighting to solve business problems, e.g., slapping a low-quality – and untested – fix onto a product and shipping it out the door to placate a customer
• Repeatedly fixing new problems yourself rather than properly training your staff to fix the problems – this is a special form known as “shifting the burden to the intervener” where you are the intervener who is inadvertently eroding the capabilities and confidence of your staff (the unintended side-effect)
• Outsourcing core business competencies rather than building internal capacity (also shifting the burden to the intervener, in this case, to the outsource provider)
• Implementing government programs that increase the recipient’s dependency on the government, e.g., welfare programs that do not attempt to simultaneously address low unemployment or low wages (also shifting the burden to the intervener, in this case, to the government)

Read more…

The simulation engine underlying STELLA® and iThink® uses numerical integration.  Numerical integration differs from the integration you may have learned in Calculus in that it uses algorithms that approximate the solution to the integration.  The two approximations currently available are known as Euler’s method and the Runge-Kutta method.  All algorithms require a finite value for DT, the integration step-size, rather than the infinitesimally small value used in Calculus.  On the surface, it may seem that the smaller DT is, the more accurate the results, but this turns out not to be true.

Compound Interest:  Euler’s Method over Runge-Kutta

To introduce Euler’s method, let’s take a look at the simple problem of compound interest.  If we have $100 that we invest at 10% (or 0.1) compounded annually, we can calculate the interest after N years by adding in the interest each year and recalculating:

1st year:  interest = $1000 × 0.1 = $100; Balance = 1000 + 100 = $1100
2nd year: interest = $1100 × 0.1 = $110; Balance = 1100 + 110 = $1210
3rd year:  interest = $1210 × 0.1 = $121; Balance = 1210 + 121 = $1331

And so on up to year N.  We have just seen the essence of how Euler’s method works.  It calculates the new change in the stock for this DT (in this case, interest) and then adds that to the previous value of the stock (Balance) to get the new value of the stock.  In this example, DT = 1 year.

By noticing we always add the existing balance in, we can instead just multiply the previous year’s balance by 1 + rate = 1 + 0.1 = 1.1:

1st year:  Balance = $1000 × 1.1 = $1100
2nd year: Balance = $1100 × 1.1 = $1210
3rd year:  Balance = $1210 × 1.1 = $1331

And so on up to year N. We can further generalize by noticing we are multiplying by 1.1 N times and thus arrive at the compound interest formula:

Balance = Initial_Balance*(1 + rate)^N

Checking this, we find our Balance at the end of year 3 is 1000*1.1^3 = $1331.  In the general case of the formula, rate is the fractional interest rate per compounding period and N is the number of compounding periods (an integer).  In our example, the compounding period is one year, so rate is the annual fractional interest rate and N is the number of years.  However, if interest is compounded quarterly (four times a year), the interest rate has to be adjusted to a per quarter rate by dividing by 4 (so rate = 0.1/4 = 0.025) and N must be expressed as the number of quarters (N = number of years*4 = 3*4 = 12 for the end of year 3).  We can use this formula in our model to test the accuracy of Euler’s method.  Note that for quarterly compounding, we would set DT = 1/4 = 0.25 years.

To explore the differences between Euler’s and Runge-Kutta, the following structure will be used for all of the examples in this post.  This structure models the compound interest problem outlined above.

The equations change for each example and can be seen in the individual model files (accessed by clicking here).  For this example, the actual value is calculated using the compound interest formula, Initial_Balance*(1 + rate)^TIME.  The approximated value is calculated by integrating rate*Approx_Balance (into Approx_Balance).

In addition to the actual and approximate values, three errors are also calculated across the model run:  the maximum absolute error, the maximum relative error, and the root-mean-squared error (RMSE).  The absolute error is:

ABS(Actual_Balance – Approx_Balance)

The relative error is:

absolute_error/ABS(Actual_Balance)

and is usually expressed as a percentage.  The RMSE is found by averaging the values of the absolute error squared, and then taking the square root of that average.

Read more…

Conveyors are useful model elements for representing pipelines or processes that take a certain amount of time to complete.  However, adding a leakage flow to a conveyor can make it difficult to initialize a model in steady-state.  The following discussion will explain how to initialize conveyors with leakage in steady-state.  Please refer to the model structure below while reading this discussion.

These additional variables will be also used:

transit_time = TRANSTIME(conveyor)
conveyor_length = transit_time/DT
leakage_fraction = the user-specified leakage fraction

Linear Leakage

The default leakage is linear in behavior.  The total amount that leaks across the length of the conveyor is directly proportional to the inflowing amount.  The leakage fraction is the constant of proportionality.  Thus, the fraction of inflowing material that makes it to the conveyor’s outflow is exactly

1 – leakage_fraction

Given the sample model structure above, to achieve equilibrium, conveyor_outflow must equal outflow.  For this to happen, we need to set the inflow as follows:

inflow = outflow/(1 – leakage_fraction)

The conveyor’s steady-state value is then:

conveyor = transit_time*inflow – (conveyor_length – 1)*leakage*DT/2

where the initial value of leakage is:

leakage = leakage_fraction*inflow

This must be calculated outside the program and entered as a constant into the conveyor as conveyors cannot be given equations (they can, however, be set to a the value of a single converter, but you must be careful how you calculate this to avoid circularity).

Exponential Leakage

Optionally, leakage can be made exponential.  The amount that leaks each DT is proportional to the amount remaining in the conveyor.  In this case, the leakage fraction is the fraction that leaks each unit of time so, for long conveyors, a lot of material can leak away.  Given the transit time, the fraction of inflowing material that makes it to the conveyor’s outflow is approximately

1 – (1 – leakage_fraction)^transit_time

Given the sample model structure above, to achieve equilibrium, conveyor_outflow must equal outflow.  For this to happen, we need to set the inflow as follows:

per_dt_no_leak = 1 – leakage_fraction*DT
inflow = outflow/(per_dt_no_leak^conveyor_length)

For steady-state, the conveyor itself must then be set as follows:

conveyor = (inflow*DT)*(1 – per_dt_no_leak^conveyor_length)/(1 – per_dt_no_leak)

I generally do not use modules to build very small models (only a couple of stocks and flows), which may then lead me to use sectors as the model grows because they are very convenient.  By the time I have three sectors, though, it starts to become clear that I should have used modules.  I will then need to convert my sector-based model into a module-based model.  Historically, I also have a number of sector-based models that are crying to be module-based.

Converting from sectors to modules is not very difficult:

1. Make sure there are no connections or flows between sectors.  Replace any of these with ghosts in the target sector.
2. In a new model, create one module for every sector.
3. Copy and paste the structure from each sector into its corresponding module.
4. Connect the modules:  At this point, the model structure has been rearranged into modules, but none of the modules are connected.  The ghosts that were in the sectors became real entities when they were pasted into the modules.  Go back to identify all of these connections and reconnect them in the module-based model.

Stepping Through a Sample Model

Let’s walk through an example.  A small sector-based model is shown below (and is available by clicking here).

This model violates what I would call good sector etiquette:  there are connectors that run between the sectors.  This is often useful in a small model such as this because it makes the feedback loops visible.  However, in a larger model, this can lead to problems such as crossed connections and difficulty in maintaining the model because sectors cannot be easily moved.

Read more…

This is the last of a three-part series on the Limits to Growth Archetype.  The first part can be accessed here and the second part here.

Last time, we explored the effects of Type 1 rivalry (rivalry between different companies in a developing market) on the Bass diffusion model by replicating the model structure.  This part will generalize this structure and add Type 2 rivalry (customers switching between brands).

Bass Diffusion with Type 1 Rivalry

To model the general case of an emerging market with multiple competitors, we can return to the original single company case and use arrays to add additional companies.  In this case, everything except Potential Customers needs to be arrayed, as shown below (and available by clicking here).

For this example, three companies will be competing for the pool of Potential Customers.  Each array has one-dimension, named Company, and that dimension has three elements, named A, B, and C, one for each company.  Although each different parameter, wom multiplier, fraction gained per $K, and marketing spend in $K, can be separately specified for each company, all three companies use the same values initially.  All three companies, however, do not enter the market at the same time.  Company A enters the market at the start of the simulation, company B enters six months later, and company C enters six months after that.

Recall that the marketing spend is the trigger for a company to start gaining customers.  Thus, the staggered market entrance can be modeled with the following equation for marketing spend in $K:

STEP(10, STARTTIME + (ARRAYIDX() – 1)*6)

The STEP function is used to start the marketing spend for each company at the desired time.  The ARRAYIDX function returns the integer index of the array element, so it will be 1 for company A, 2 for company B, and 3 for company C.  Thus, the offsets from the start of the simulation for the launch of each company’s marketing campaign are 0, 6, and 12, respectively.

This leads to the following behavior:

Note that under these circumstances, the first company to enter the market retains a leadership position.  However, companies B and C could anticipate this and market more strongly.  What if company B spent 50% more and company C spent 100% more than company A on marketing that is similarly effective?  This could be modeling by once again changing the equation for marketing spend in $K, this time to:

STEP(10 + (ARRAYIDX() – 1)*5, STARTTIME + (ARRAYIDX() – 1)*6)

Read more…

This is part two of a three part series on Limits to Growth.  Part one can be accessed here and part three can be accessed here.

In part one of this series, I explained the Limits to Growth archetype and gave examples in epidemiology and ecology. This part introduces the Bass diffusion model, an effective way to implement the capture of customers in a developing market. This is also used to implement what Kim Warren calls Type 1 rivalry in his book Strategy Management Dynamics, that is, rivalry between multiple companies in an emerging market.

The Bass Diffusion Model

The Bass diffusion model is very similar to the SIR model shown in part one. Since we do not usually track customers who have “recovered” from using our product, the model only has two stocks, corresponding loosely to the Susceptible and Infected stocks. New customers are acquired through contact with existing customers, just as an infection spreads, but in this context this is called word of mouth (wom). This is, however, not sufficient to spread the news of a good product, so the Bass diffusion model also includes a constant rate of customer acquisition through advertising. This is shown below (and can be downloaded by clicking here).

The feedback loops B1 and R are the same as the balancing and reinforcing loops between Susceptible and Infected in the SIR model. Instead of an infection rate, there is a wom multiplier which is the product of the Bass diffusion model’s contact rate and the adoption rate. If you are examining policies related to these variables, it would be important to separate them out in the model.

The additional feedback loop, B2, starts the ball rolling and helps a steady stream of customers come in the door. If you examine the SIR model closely, you will see that the initial value of Infected is one. If no one is infected, the disease cannot spread. Likewise, if no one is a customer, there is no one to tell others how great the product is so they want to become customers also. By advertising, awareness of the product is created in the market and some people will become customers without having encountered other customers who are happy with the product.

The behavior of this model is shown below. Note it is not different in character from the SIR model or the simple population model.

Read more…

This is the first of a three-part series on the Limits to Growth Archetype.  The second part can be accessed here and the third part here.

The Limits to Growth Systems Archetype, also known as Limits to Success, combines growth with an exogenous or endogenous limit.  This Systems Archetype was formally identified in Appendix 2 of The Fifth Discipline by Peter Senge (1990), but made its first prominent appearance in World Dynamics by Jay Forrester (1971) and then The Limits to Growth by Meadows, Meadows, Randers, and Behrens (1972).  The Causal Loop Diagram (CLD) is shown below.

Real growth processes have inherent limits to growth.  Identifying these limits can help avoid problems in the future, whether the problem is overpopulation, increasing demand for a product that cannot be met, or growing a business in a mature market.  When growth is desired, but limited, it is always better to find ways to increase the limit before pushing for more growth.  Excessive growth in the face of a limit often leads to collapse.  Driving the system to the point of collapse can erode the ability to continue after the collapse, for example, by reducing the production capability of a piece of farmland or destroying the reputation of a company.

Classic examples of limits to growth include:

• The collapse of the deer population on the Kaibab plateau and on St. Matthew Island due to overpopulation and the attendant overgrazing of their habitat
• The overshoot and collapse of the human population on Easter Island
• Overgrazing in the Sahel region of Africa by cattle herders
• Overfishing of the oceans by fishermen
• The collapse of People Express due to sharp customer growth combined with slow personnel growth
• The sharp exodus of America Online subscribers after an intense marketing campaign increased the number of subscribers far beyond their capacity
• The contraction of the world economy in 2008 due to limiting oil supplies
• The productivity of staff deteriorating as a company grows, due to increased interactions and reporting overhead
• Business growth limited by the size of the potential market
• Yeast cells in the fermentation process, who suffer from both the loss of exogenously supplied sugar and the increase of endogenously produced pollution

Read more…