In this blog we will cover Discrete-Time Markov chains (or simply Markov chains) in R. We assume you have a basic understanding of following tasks :
Although we will be revisiting these concepts but having a basic understanding will help you get the most out of this blog.
Introduction
Let's first define a Stochastic process:
A Stochastic process, Xt is a set of ordered random variables defined for each time t in some set J (time-domain). The time-domain may be discrete or continuous. The set of values that the Xt is capable of taking is called the state space S. The state-space may also be discrete or continuous.
For example, we might model the number of wickets fallen after every over in a cricket match (20 overs) by a stochastic process. Here both the time-domain and the state space will be discrete. J = {0, 1, 2,....., 20) overs and S = {0, 1, 2,....., 10} wickets.
A Markov chain is a stochastic process in discrete-time and with a discrete state space satisfying the Markov property. The Markov property tells us that if we know the current value of the process, we have all the information required to determine the probabilities for future values. Knowing the historical values of the process as well would not make any difference.
An example is a no claims discount (NCD) model run by motor insurance companies. A company might offer discounts of 0%, 25%, 40%, and 60% of the full premium.
Let's consider some examples other from Insurance.
Google's famous PageRank algorithm, typing word prediction on keypads, weather prediction, and many more.
One of the most famous use cases of Markov Chains is Google’s famous PageRank algorithm. Google has recorded all the websites and created a transition matrix on a very large scale. Stages are represented by the sites the users are currently on, and the transition matrix represents the user's probability of moving from one site to another site.
One way of representing a Markov chain is by its transition graph, and another way is by its transition probability matrix. For the above example of NCD model, the transition probability matrix will be as follows:
where pij denotes the transition probability from state i to state j, where i,j = 1, 2, 3, 4.
The long-term distribution of a Markov chain
πj, j ∈ S is a stationary probability distribution for a Markov chain with transition probability matrix P if the following conditions hold for all j in S:
πj = Σ πipij for i ∈ S
πj ≥ 0
Σ πj = 1 for j ∈ S
In a compact form π = πP.
This means that if the chain reaches the distribution π at some time n, then the distribution of Xt will be π for all subsequent t ≥ n. That is, the statistical properties of the process do not change over time, so the chain is called a stationary process.
Results for Stationary distribution
A Markov chain with a finite state space has at least one stationary probability distribution.
An irreducible Markov chain with a finite state space has a unique stationary probability distribution.
Irreducibility - A Markov chain is said to be irreducible if any state j can be reached from any other state i.
An irreducible, aperiodic Markov chain with a finite state space will settle down to its unique stationary distribution in the long run.
Periodicity - A state i is said to be periodic with period d > 1 if a return to i is only possible in a number of steps that is a multiple of d. A state is said to be aperiodic if it is not periodic.
If a Markov chain is irreducible then all its states have the same period (or all are aperiodic).
A Markov chain is said to be aperiodic if all its states are aperiodic.
A simple No Claims Discount (NCD) policy
In the No Claims Discount (NCD) system in motor insurance, the premium charged depends on the driver's claim record. This is the main example of Markov chains.
Let's first consider a simple No Claims Discount (NCD) model.
A motor insurance company offers the following discounts to motorists on their annual premium:
Level 1 0% discount
Level 2 25% discount
Level 3 40% discount
After a claim‐free year, probability 0.75, policyholders move up a level or stay at Level 3.
After a year with one or more claims, probability 0.25, policyholders move down a level or stay at Level 1.
Transition graph:
Transition matrix:
Is this a Markov process? The process has the Markov property since the probability of moving on to the next level only depends on the current level (so it is not dependent on the past history of the process). So the process is a Markov process.
State-space:
S = {Level1, Level2, Level3}
Discrete state-space
Time-domain:
J = {0, 1, 2,......} years
Discrete time-domain
Is this a Markov chain? Since the process has both discrete state-space and discrete time-domain and the process also satisfies the Markov property. Hence this is a Markov Chain.
Is the chain irreducible? Since any state j can be reached from any other state i (not necessarily in one step), the chain is irreducible, where i, j ∈ S.
Is the chain aperiodic? Since the chain is irreducible, all its states have the same period. Consider state 1 (Level 1), since state 1 has an arrow back to itself, it is possible to return to state 1 in 1 steps. Thus, it is aperiodic and hence the chain is also aperiodic.
Using the results of the stationary distribution, the chain will settle down to its unique stationary distribution in the long run.
Setting up the one‐step transition probability matrix for this process in R:
We will use the matrix function to set this matrix up in R. One efficient way is to first create a 3 X 3 matrix full of zeros and then later populate it:
P = matrix(0, nrow = 3, ncol = 3)and then, populating the matrix with non-zero probabilities:
P[1,1] = 0.25
P[1,2] = 0.75
P[1,3] = 0
P[2,1] = 0.25
P[2,2] = 0
P[2,3] = 0.75
P[3,1] = 0
P[3,2] = 0.25
P[3,3] = 0.75Displaying the matrix:
P
[,1] [,2] [,3]
[1,] 0.25 0.75 0.00
[2,] 0.25 0.00 0.75
[3,] 0.00 0.25 0.75We can use rowSums function to check that the probabilities in each row of the matrix add up to 1. (In this case it is straightforward that the probabilities in each row add up to 1 but this is a useful function when dealing with larger matrices).
rowSums(P)
[1] 1 1 1Consider the following situation: A new policyholder is currently on the 0% discount level (Level 1). The motor insurance company needs to calculate the probabilities for the discount rate this policyholder will be receiving after 5 years.
This can be calculated as follows:
We first start by creating an initial state vector for a new policyholder.
initial_state = c(1, 0, 0)Then calculating the transition probability matrix after 5 years, P^5 by successively multiplying P five times:
P5 = P%*%P%*%P%*%P%*%PNote: Here we are using %*% for matrix multiplication instead of ^ symbol which simply raises each element of the matrix to some power.
Finally, to calculate the distribution after 5 years we simply multiply the initial_state vector with P^5 (transition probability matrix after 5 years).
initial_state %*% P5
[,1] [,2] [,3]
[1,] 0.083 0.249 0.668Thus, the probabilities for the new policyholder’s discount rate after 5 years are given as follows:
Level 1 - 0%: 8.3%
Level 2 - 25%: 24.9%
Level 3 - 40%: 66.8%
Alternatively, we can code a function that will speed up things for us because we often need an n‐step transition matrix.
Pn = function(P, n) {
M = diag(nrow(P))
for (i in 1:n) {
M = M %*% P
}
M
}The function work as follows:
The function takes two arguments - the probability matrix P and the argument n to specify the n-steps.
nrow(P) computes the number of rows of the probability matrix P.
The diag() function then gives us the identity matrix of the same dimension as that of probability matrix P.
Then we ran a for loop from 1 to n (n-steps) and multiply the diagonal matrix with the probability matrix and update it for each loop.
We can then use this function to calculate P^5:
P5 = Pn(P,5)
P5
[,1] [,2] [,3]
[1,] 0.083 0.249 0.668
[2,] 0.083 0.223 0.694
[3,] 0.074 0.232 0.694
Finally, multiplying the initial_state vector with P5 gives us the distribution after 5 years.
initial_state %*% P5
[,1] [,2] [,3]
[1,] 0.083 0.249 0.668Same as above!
Let us say the motor insurance company is now interested in knowing the distribution for the discount rate the new policyholder will be receiving after 20 years.
For this, we simply need to compute initial_state %*% P^20 and we use the function created above to calculate P^20.
P20 = Pn(P,20)
P20
[,1] [,2] [,3]
[1,] 0.07692 0.2308 0.6923
[2,] 0.07692 0.2308 0.6923
[3,] 0.07692 0.2308 0.6923 We can see that the probabilities in each row of P20 are very similar. This indicates that the process converges to a stationary distribution in the long‐run with these probabilities.
initial_state %*% P20
[,1] [,2] [,3]
[1,] 0.07692 0.2308 0.6923Thus, the probabilities for the new policyholder’s discount rate after 20 years are given as follows:
Level 1 - 0%: 7.692%
Level 2 - 25%: 23.077%
Level 3 - 40%: 69.231%
The company is also interested in knowing the expected number of motorists at each level after 2 years from a group of 100 motorists, all initially at Level 1.
We now define the initial state as follows:
initial_state = c(100, 0, 0)The two-step transition probability matrix can be calculated using the above-defined function Pn or by P %*% P.
P2 = Pn(P, 2)
P2
[,1] [,2] [,3]
[1,] 0.2500 0.1875 0.5625
[2,] 0.0625 0.3750 0.5625
[3,] 0.0625 0.1875 0.7500Thus, the expected number of motorists at each level after 2 years are:
initial_state %*% P2
[,1] [,2] [,3]
[1,] 25 18.75 56.25That is, it is expected that the number of motorists in Level 1, 2, and 3 after 2 years will be approximately 25, 19, and 56 respectively.
Now, consider the initial premium of $100 for the motorists. The motorists at Level 2 will get a 25% premium discount and those in Level 3 will get a 40% premium discount. Thus, a motorist at Level 1 will pay the full premium of $100 per annum, a motorist at Level 2 will pay a premium of $75 per annum, and a motorist at Level 3 will pay a premium of $60 per annum.
The company is now interested in calculating its expected premium income over the next 3 years from a new motorist.
We will write a for-loop that adds up the premium income over the next n years for a given starting distribution. Before that, we create a row vector of discount levels, the current position of the new motorist and set the premium income to be 0:
position = c(1,0,0)
discount = c(0,0.25,0.40)
prem.inc = 0We start a loop for 3 years:
for (j in 1:3) {
prem.inc = prem.inc + 100*sum(position*(1-discount))
position = position %*% P
}The loop work as follows:
For each j in the loop, we add that year’s premium income to the total
We did this by multiplying the position of the motorist with (1 - discount) level and summing it and multiplying with the premium of $100,
and then update the position of drivers.
Finally printing the final premium income after 3 years for a new motorist:
prem.inc
[1] 254.0625What R has done here is to take the first year’s premium of $100 and then apply P to the new driver:
So, in the second year, with 75% of drivers receiving a 25% discount, the premium income is:
= 0.25 X 100 + 0.75 X 75 = 81.25
So, R adds this to the $100, giving us a total so far of $181.25.
We apply P again:
So, in the third year, with 18.75% of drivers receiving a 25% discount, and 56.25% of drivers receiving a 40% discount the premium income is:
= 0.25 X 100 + 0.1875 X 75 + 0.5625 X 60 = 72.8125
Finally, R adds this to the $181.25 to get $254.0625.
Using for loop we got the same premium income of $254.0625.
Thus, the company is expected to earn $254.0625 premium in 3 years from a new motorist.
That's all for this blog but not for the topic! In the next blog, we will be exploring markovchain package in R and see how it can be used for making things easier for us and we will also see its uses in modelling and simulating Markov chains.
Thank you for reading the post. We hope you enjoyed it! Do let us know in the comments if you enjoyed reading this post. We would be very grateful if you’d help it spread by sharing it with your friends and network. Thank you!
—Abhishek Goel and Bhavyadeep Gupta
We both are currently actuarial students and the above blog is to the best of our knowledge. We do not claim to be the experts in the field. We are just sharing our learning with others and this is an effort to help other students learn R. Your feedback and suggestions are welcomed.
تعليقات