When I started working on this How-To on building a simple Markov chain (a useful component of model-building), I came across this great visualization at Setosa Blog. Clearly, after this and the really incredible Segregation visualization that Iza posted about, I’m starting to feel the need to step up my data viz game. So before you read on, have a look at the visualization because it does a really good job of explaining the fundamentals and gives you a chance to play with a Markov chain yourself. I’ll just briefly cover what Markov chains are, and then we’ll get into using Markov chains in an ABM using NetLogo.

**What is a Markov chain?**

A Markov chain is a way to model how a system changes from one state to another over time. Imagine a system with two states: A and B. If the system is in state A, there is a probability that over the next time step, the system will transition to state B (with the inverse probability that the system will remain in state A). There are also probabilities that a system in state B will transition to state A. Since the system we’re modeling has to exist in one of the two states, the transition probabilities for each state should add up to 1.00. These are sometimes visualized using a simple network graph, like so:

Markov chains are useful for dealing with phenomena that are auto-correlated: that is, the future state of the phenomenon depends on its current state. When this is true, we may not be able to accurately model the phenomenon as a random probability draw. To use this method, we need to know what the current state of the phenomenon is, and what the likelihood is of the system transitioning to another state. The weather is often used as an example. In ecological modeling, we often want to model rainfall because pretty much every living thing on earth depends to some degree on the frequency and regularity of rainfall. Let’s look at the rainfall data for Boston, Massachusetts in April of last year (courtesy of Weather Underground).

It rained 16 out of 30 days in April, or 53.33% of the time. We could simulate April rain frequency in Boston by simply performing 30 checks against a probability of 53.33%. If you wanted to do that in R, we could do the following:

`prob<-0.5333 `

`checks<-(runif(30, 0, 1)) <= prob`

Here, we have two variables: one expressing the transition probability, and the other running the routine of generating thirty random numbers between 0 and 1 and checking whether they are less than or equal to that probability. If we enter **checks** into the command prompt, it would give us the following output:

`TRUE FALSE TRUE TRUE TRUE FALSE TRUE TRUE TRUE FALSE FALSE TRUE FALSE TRUE FALSE TRUE FALSE FALSE TRUE TRUE TRUE TRUE TRUE FALSE TRUE FALSE TRUE FALSE TRUE FALSE`

In this simulated dataset, TRUE is a rainy day and FALSE is a clear day. That looks OK, but does rain really work that way? Do we really get long spans where it rains on alternating days? One thing we notice from the Boston rainfall data is that rainy days don’t seem to occur sporadically; they occur in groups in between spans of good weather. This probably makes intuitive sense; if it starts raining today, there’s a better-than-average chance it might rain tomorrow. That’s because the process of rain clouds passing over our heads is auto-correlated. However, we can use this data to estimate how likely it is for the weather to go from dry to wet (DTW), or vice versa (WTD). We calculate these values using the following equations:

DTW = number of dry days followed by a wet day / total number of dry days = 5 / 14 = 35.71%

WTD = number of wet days followed by a dry day / total number of wet days = 5 / 16 = 31.25%

More data would be better, and Weather Underground has Boston weather going back about 95 years. But with these transition probabilities, we could build a Markov chain that simulates the sequence of rainy versus clear days in April in Boston in a way that captures some of that auto-correlation.

**Making it rain in NetLogo**

Let’s say we wanted to model grass growth in NetLogo, using a Markov chain like the one above to simulate daily rainfall. First, we’ll open a new model, head over to the **Code** tab, and create a global variable for the **weather. **We can also consider grass as a component of the patches in the model, and the grass could be in various growth stages, which we’ll call **growth**.

Next, we’ll set the model up so that we start off in one of two weather conditions: “**WET”** or “**DRY”**. We’ll also set up the grass so that growth is randomly distributed between 0 and 10 among the patches (we code this as **random 11** to account for 0), and we’ll color the patches according to their growth using the **scale-color** command (darker = more growth, lighter = less growth).

Then, in the **go** procedure, we’ll have two functions: **markov-rain** and **grow-grass**, that will repeat during the model at each time step, or **tick**, until a limit of 500 ticks is reached.

The **markov-rain** procedure will determine the transitions between **“WET”** and **“DRY”** weather, and we’ll use variable placeholders **WTD** for the wet-to-dry transition, and **DTW** for the dry-to-wet transition. This might throw up an error about these variables not existing, but we’ll get to that later.

What’s happening here? To use a Markov chain in NetLogo, we need to know what the current state of the process is, in this case the current **weather**, and then apply our transition probabilities to determine whether the state has changed. Here, we use the **ifelse** command to determine whether it is currently **“WET”** or **“DRY”**, and then determine the next state using the appropriate transition probability (**WTD** if conditions are currently **“WET”**, **DTW** if conditions are currently **“DRY”**). If a random number draw between 0 and 1.000 is less than whatever the value for the transition is, then the **weather** will change to its opposite state; otherwise, it will remain the same. The **grow-grass** procedure will again use **ifelse** to determine the **weather**. If it is **“WET”**, any grass which has not reached the maximum **growth** level of 10 will grow by one unit. If it is **“DRY”**, grass which has not reached the minimum **growth** level of 0 will grow by one unit.

The last step is to switch over to the Interface tab and create sliders for the **WTD** and **DTW** variables. I set these up from 0 to 1.0 at an interval of 0.1. I also added the **setup** and **go** buttons as well, making sure to tick the **Forever?** option for the **go** button so that it will run repeatedly until the time limit is reached. Finally, I added a plot which records the **mean [ growth ] of patches**. The interface should look something like this:

**Running the model**

When both the **WTD** and **DTW** are set to 50%, the rainfall effectively follows a random walk between dry and wet, as we might expect. This produces grass growth that varies randomly over time.

If we use the sliders to change the conditions so that it is likelier to transition from dry to wet (**DTW** = 80%) than wet to dry (**WTD**= 20%), the outcome is fairly predictable: wetter conditions on the whole, meaning greater vegetation growth (top image). When the transition is the other way around (**DTW** = 20%, **WTD** = 80%), conditions are drier which means less grass growth (bottom image).

But what about when both transitions are less likely to occur? Let’s set them both at 30%, not far off from our Boston scenario:

Under this setting we get dramatic sways between dry and wet conditions, causing rapid growth and decline of vegetation. What about if both transitions are higher? Let’s set both transitions to 80%:

When both transitions are set to 80%, we get fewer dramatic sways in vegetation growth. Why does this happen? The answer lies in the likelihood of transitioning and the buffer provided by the growth factor. In the 30%-30% scenario, the likelihood of transitioning overall is low, but when it does happen, it’s more likely to stay in whatever state it transitions to, swinging the vegetation to one extreme or the other fairly rapidly. Under the 80%-80% scenario, there is a greater likelihood that a wet day will be quickly offset by a dry day (and vice versa), which has a balancing effect over the short-term. This is an interesting behavior that may not be intuitive without building and running the model first.

**Going further**

Of course, this is a very simple model, and I sincerely hope no one would actually attempt to predict the weather or vegetation growth with it (DISCLAIMER: No, really, don’t do it). How could it be improved? Well, for starters, while we model rain/not rain in binary terms here, we know that the intensity of rain can vary tremendously. Since there is no limit to the number of states we can use, we could add further numbers of states and transitions to model the likelihood of going from drizzle to deluge and all points in between. Furthermore, the transition probabilities are likely to change over the course of the year, relative to seasonal changes. There are also likely to be changes from year to year or decade to decade relative to global climatic cycles. Adding these different levels of change could be accomplished by treating each smaller scale time period as nodes within a larger Markov chain with its own transition probabilities.

But we might also use Markov chains for modeling adaptive agent behaviors as well. Oftentimes, we want agents that don’t follow the same script by rote each time step. For instance, an agent might eat at its favorite restaurant x% of the time, but y% of the time it might try something new. How do different settings of x and y affect the ranging patterns of the agent? What if the agent changes the likelihood of eating at particular restaurants based on dining experiences? What determines how likely a restaurant is to be visited? Each agent could develop an entire decision-making schema and subsequent evolving ranging patterns through the process of exploring and adapting a simple Markov chain like the one presented here.

**More on Markov chains**

One of the interesting features of a Markov process with a finite number of states and fixed transition probabilities (like our rainfall model) is that, over time, it converges on a distribution of outcomes. This can be really useful for predicting the future state of the system, but also for determining where changing the transition probabilities is likely to have the greatest effect. There’s a lot more to learn about Markov chains and what kinds of processes are best represented by them. Here are a few things to get you started:

- The historical background on Markov chains gets explained like a pro in this video from Khan Academy.
- This set of short videos from Scott Page’s Model Thinking course does a good job of providing an overview, and explains the Markov Convergence Theorem in more detail.
- This freely available book has a good overview of how Markov chains work mathematically in Chapter 11.
- This paper by Izquierdo and colleagues gives a very thorough treatment of computer simulation and Markov chain analysis.

*Featured image: A Markov chain of how my PhD thesis is coming along…*

UPDATE 18 Feb 2015: Per request, the code for the model can be downloaded here: markov model

UPDATE 10 Aug 2015: Full model available here