Agenda

• Packages for Markov Chain
• Transition Probability Matrix
• Creating Markov Chain R object (CoCo Example)
  • From TPM
  • From State Sequence
  • Absorbing Markov Chain
• Markov Chain Operations
  • TPM after n steps
  • State Probability after n steps
  • Steady State Probability

Agenda

• Absorbing Markov Chains (Airwaves India Example)
  • Fundamental Matrix
  • Time to Absorption
  • Absorption Probability

Packages for Markov Chain

library(expm) #Matrix Exponential, Log, etc
library(diagram) #Functions for visualising simple graphs 
library(markovchain)

Calculating TPM from State Sequences

rain dataset from markovchain package
- Probability of 0cm, 1-5cm, 6+ cm of rain on the island of Alofi

data(rain) 
head(rain)

##   V1 rain
## 1  3   6+
## 2  2  1-5
## 3  2  1-5
## 4  2  1-5
## 5  2  1-5
## 6  2  1-5

Calculating TPM from State Sequences

rain dataset from markovchain package

rain.matrix <- createSequenceMatrix(rain$rain)

rain.matrix

##       0 1-5  6+
## 0   362 126  60
## 1-5 136  90  68
## 6+   50  79 124

Understand the sequence data

sequence <- c("a", "b", "a", "a", "a", "a", "b", "a", "b", "a", "b", "a", "a", 
              "b", "b", "b", "a")
sequenceMatr <- createSequenceMatrix(sequence)
sequenceMatr

##   a b
## a 4 5
## b 5 2

#verifyMarkovProperty(sequence)

Calculating TPM from State Sequences

rain dataset from markovchain package

prop.table expresses table Entries as Fraction of Marginal Table

rain.tpm <- prop.table(rain.matrix, margin = 1)

show(rain.tpm)

##             0       1-5        6+
## 0   0.6605839 0.2299270 0.1094891
## 1-5 0.4625850 0.3061224 0.2312925
## 6+  0.1976285 0.3122530 0.4901186

Markov Chain Object from State Sequence

markovchainFit(): Given a sequence of states arising from a stationary state, it fits the underlying Markov chain distribution using either MLE or other methods.

Use fitHigherOrder() to fit a higher order markov chain

rain.mc <- markovchainFit(data = rain$rain)

When the “mle”, “bootstrap” or “map” method is used, the lower and upper confidence bounds are returned along with the standard error: - Understand MLE from this link - For explanation on Markov Chain logLikelihood statistics refer this

rain.mc

## $estimate
## MLE Fit 
##  A  3 - dimensional discrete Markov Chain defined by the following states: 
##  0, 1-5, 6+ 
##  The transition matrix  (by rows)  is defined as follows: 
##             0       1-5        6+
## 0   0.6605839 0.2299270 0.1094891
## 1-5 0.4625850 0.3061224 0.2312925
## 6+  0.1976285 0.3122530 0.4901186
## 
## 
## $standardError
##              0        1-5         6+
## 0   0.03471952 0.02048353 0.01413498
## 1-5 0.03966634 0.03226814 0.02804834
## 6+  0.02794888 0.03513120 0.04401395
## 
## $confidenceLevel
## [1] 0.95
## 
## $lowerEndpointMatrix
##             0       1-5         6+
## 0   0.6034754 0.1962346 0.08623909
## 1-5 0.3973397 0.2530461 0.18515711
## 6+  0.1516566 0.2544673 0.41772208
## 
## $upperEndpointMatrix
##             0       1-5        6+
## 0   0.7176925 0.2636194 0.1327390
## 1-5 0.5278304 0.3591988 0.2774279
## 6+  0.2436003 0.3700387 0.5625151
## 
## $logLikelihood
## [1] -1040.419

Markov Chain Object from State Sequence

summary(rain.mc$estimate)

## MLE Fit  Markov chain that is composed by: 
## Closed classes: 
## 0 1-5 6+ 
## Recurrent classes: 
## {0,1-5,6+}
## Transient classes: 
## NONE 
## The Markov chain is irreducible 
## The absorbing states are: NONE

Markov Chain Object from State Sequence

rain.mc$estimate@transitionMatrix

##             0       1-5        6+
## 0   0.6605839 0.2299270 0.1094891
## 1-5 0.4625850 0.3061224 0.2312925
## 6+  0.1976285 0.3122530 0.4901186

states(rain.mc$estimate)

## [1] "0"   "1-5" "6+"

Plotting the Markov Chain

plot(rain.mc$estimate, package = 'diagram')

#plotmat(rain.mc)

Calculating Transition Matrix from Raw Data

The actual data may look like this:

dat = data.frame(replicate(250,sample(c("A", "B", "C","D"), size = 6, replace=TRUE)))

head(dat)

##   X1 X2 X3 X4 X5 X6 X7 X8 X9 X10 X11 X12 X13 X14 X15 X16 X17 X18 X19 X20
## 1  B  A  B  C  C  B  C  C  C   B   B   D   B   D   C   B   A   C   A   D
## 2  D  C  C  A  A  B  D  D  A   B   A   A   D   C   B   B   B   D   B   C
## 3  A  C  D  C  B  A  A  B  B   C   B   B   C   B   D   B   C   C   B   B
## 4  C  A  D  D  B  D  A  D  C   B   D   C   B   C   D   A   A   B   C   A
## 5  D  C  A  A  D  B  B  A  B   B   D   D   A   D   C   C   A   B   D   D
## 6  C  D  C  B  C  A  B  D  B   B   D   C   D   A   D   A   D   C   A   D
##   X21 X22 X23 X24 X25 X26 X27 X28 X29 X30 X31 X32 X33 X34 X35 X36 X37 X38
## 1   B   D   A   C   B   D   A   A   C   A   D   B   D   B   A   D   D   D
## 2   D   A   B   B   D   C   B   A   C   C   D   B   B   A   A   C   B   C
## 3   C   B   A   D   B   D   B   A   B   D   A   C   B   A   C   C   B   B
## 4   B   C   B   D   C   C   D   B   B   D   C   C   C   D   B   D   C   B
## 5   B   D   A   B   A   D   C   B   A   D   C   A   D   D   A   B   B   C
## 6   B   A   C   A   B   A   B   B   C   B   C   C   A   A   B   A   B   D
##   X39 X40 X41 X42 X43 X44 X45 X46 X47 X48 X49 X50 X51 X52 X53 X54 X55 X56
## 1   B   C   B   C   B   A   D   D   D   C   B   A   C   A   D   D   D   A
## 2   C   A   D   C   A   D   D   D   D   B   A   C   B   D   C   D   A   B
## 3   A   A   B   D   A   A   D   A   B   A   C   B   B   A   A   A   D   D
## 4   B   D   C   C   B   D   D   D   A   B   B   C   A   C   C   B   C   C
## 5   A   C   D   C   A   D   B   C   A   B   C   C   D   C   A   C   A   C
## 6   D   C   D   A   C   D   A   D   C   B   D   B   B   A   B   D   B   B
##   X57 X58 X59 X60 X61 X62 X63 X64 X65 X66 X67 X68 X69 X70 X71 X72 X73 X74
## 1   C   A   B   A   D   C   B   A   B   A   B   A   D   D   A   C   B   A
## 2   B   A   A   C   B   B   A   C   D   B   A   C   C   C   C   C   B   B
## 3   D   A   A   B   B   C   C   D   D   D   C   C   D   D   D   C   D   B
## 4   A   A   B   A   B   B   C   B   D   D   A   C   D   B   A   D   C   C
## 5   B   D   A   C   B   A   C   B   D   B   B   D   B   D   B   D   C   A
## 6   B   B   C   C   C   C   C   D   D   C   B   D   B   C   D   B   D   B
##   X75 X76 X77 X78 X79 X80 X81 X82 X83 X84 X85 X86 X87 X88 X89 X90 X91 X92
## 1   A   A   B   C   D   B   B   C   B   C   B   D   D   A   C   A   A   B
## 2   A   C   B   C   D   B   C   C   A   C   B   D   B   D   C   A   D   C
## 3   A   A   D   C   B   D   C   D   C   C   A   B   B   A   B   A   B   C
## 4   B   D   D   D   C   D   B   D   C   A   C   B   A   A   B   B   D   D
## 5   A   B   D   D   D   A   D   B   D   B   A   A   D   D   D   D   C   C
## 6   C   B   C   A   B   B   A   A   C   A   A   A   D   A   C   C   B   C
##   X93 X94 X95 X96 X97 X98 X99 X100 X101 X102 X103 X104 X105 X106 X107 X108
## 1   C   B   A   B   C   B   D    C    B    A    A    A    D    C    D    A
## 2   B   D   B   A   D   A   C    D    A    B    A    D    B    C    B    A
## 3   B   C   C   B   B   A   D    A    D    A    B    B    A    D    B    C
## 4   C   C   A   A   B   C   D    D    A    A    C    C    D    D    B    A
## 5   C   C   B   A   B   C   A    C    B    D    B    A    A    A    C    B
## 6   C   A   B   C   B   B   D    D    D    A    B    C    A    C    C    C
##   X109 X110 X111 X112 X113 X114 X115 X116 X117 X118 X119 X120 X121 X122
## 1    D    A    B    A    D    B    A    C    C    C    C    A    D    C
## 2    B    C    B    D    C    C    A    D    C    C    D    D    A    D
## 3    A    D    D    C    A    D    C    B    D    D    B    A    A    D
## 4    A    D    A    D    D    D    D    B    C    A    B    C    C    B
## 5    D    C    A    A    D    C    B    A    C    B    C    C    A    A
## 6    D    D    D    D    C    A    D    D    A    D    C    B    C    A
##   X123 X124 X125 X126 X127 X128 X129 X130 X131 X132 X133 X134 X135 X136
## 1    D    D    B    D    C    D    D    D    D    A    B    C    D    D
## 2    C    C    D    D    A    C    B    B    B    B    A    D    C    A
## 3    D    B    A    A    B    C    C    D    A    A    B    C    A    D
## 4    A    C    C    B    B    B    B    C    B    D    A    B    C    C
## 5    B    A    D    D    C    A    C    D    C    A    D    A    D    A
## 6    C    A    D    A    C    C    C    A    A    A    B    C    A    C
##   X137 X138 X139 X140 X141 X142 X143 X144 X145 X146 X147 X148 X149 X150
## 1    D    B    A    D    A    C    C    A    D    D    A    D    C    B
## 2    B    A    C    C    A    B    B    A    D    C    D    C    B    D
## 3    D    C    A    C    C    A    D    B    C    C    A    B    C    B
## 4    B    C    A    A    B    C    D    B    D    D    C    B    D    D
## 5    C    C    C    D    C    C    C    B    C    C    C    B    B    B
## 6    C    D    B    A    C    D    D    B    D    A    C    B    B    D
##   X151 X152 X153 X154 X155 X156 X157 X158 X159 X160 X161 X162 X163 X164
## 1    B    D    C    C    D    B    D    C    A    B    A    A    D    D
## 2    A    D    A    B    D    C    C    C    A    B    D    A    C    D
## 3    A    C    B    D    D    B    B    D    C    A    A    B    C    A
## 4    A    B    D    D    B    A    C    A    B    A    C    C    B    B
## 5    A    C    C    A    C    C    D    C    B    D    A    A    D    A
## 6    C    A    D    C    B    A    C    C    A    D    B    C    D    A
##   X165 X166 X167 X168 X169 X170 X171 X172 X173 X174 X175 X176 X177 X178
## 1    A    A    C    A    D    A    B    C    C    C    D    B    D    A
## 2    B    D    D    B    A    C    B    A    C    D    D    D    D    A
## 3    B    B    D    A    C    C    A    B    C    A    D    D    B    A
## 4    C    A    A    B    D    C    C    B    D    A    D    B    B    C
## 5    D    A    D    C    C    A    A    D    C    C    C    B    C    C
## 6    B    A    B    D    A    C    D    D    B    A    A    A    B    C
##   X179 X180 X181 X182 X183 X184 X185 X186 X187 X188 X189 X190 X191 X192
## 1    A    D    B    B    C    C    B    B    D    C    A    C    A    B
## 2    C    D    D    A    B    C    C    B    A    C    D    D    A    B
## 3    D    A    D    C    B    C    A    D    B    A    B    B    D    D
## 4    C    C    A    D    B    D    D    B    B    D    C    D    A    B
## 5    C    C    B    C    B    A    B    A    A    B    B    A    C    C
## 6    C    B    D    D    A    B    A    C    A    D    C    A    C    B
##   X193 X194 X195 X196 X197 X198 X199 X200 X201 X202 X203 X204 X205 X206
## 1    C    C    A    C    B    C    C    A    A    A    D    D    B    B
## 2    D    C    D    C    A    D    B    D    C    A    D    D    D    C
## 3    D    A    D    C    D    C    D    B    C    A    A    C    C    B
## 4    D    A    D    B    B    C    A    A    B    C    B    B    C    C
## 5    C    C    D    C    C    B    A    D    C    C    D    C    D    A
## 6    D    D    B    A    A    B    D    D    C    C    D    C    B    B
##   X207 X208 X209 X210 X211 X212 X213 X214 X215 X216 X217 X218 X219 X220
## 1    D    A    A    C    D    D    C    A    D    D    B    B    D    A
## 2    D    B    D    D    D    C    D    B    D    C    D    C    B    B
## 3    B    A    B    B    D    A    B    B    B    B    C    B    D    D
## 4    B    B    A    A    C    B    C    A    C    A    B    C    A    C
## 5    C    D    C    A    D    A    A    C    A    A    D    B    C    A
## 6    C    A    D    A    D    A    B    D    D    C    B    C    A    D
##   X221 X222 X223 X224 X225 X226 X227 X228 X229 X230 X231 X232 X233 X234
## 1    C    B    C    B    B    C    A    B    D    A    A    C    D    B
## 2    D    A    B    C    A    D    C    B    B    A    C    B    B    B
## 3    D    D    A    A    A    C    C    B    B    C    C    B    C    D
## 4    A    A    A    B    D    D    C    C    B    C    C    B    C    A
## 5    D    C    D    A    A    D    D    B    C    A    C    B    B    B
## 6    C    C    D    C    B    B    C    C    D    D    B    C    C    C
##   X235 X236 X237 X238 X239 X240 X241 X242 X243 X244 X245 X246 X247 X248
## 1    C    B    D    D    A    B    A    A    B    C    A    A    D    C
## 2    C    D    D    C    C    A    D    A    D    C    D    C    B    A
## 3    D    A    C    B    D    C    B    A    D    D    D    B    C    C
## 4    B    D    D    B    B    D    A    A    A    B    B    A    D    A
## 5    A    A    C    A    B    D    D    C    B    B    C    C    D    C
## 6    A    C    B    D    C    C    B    D    C    C    A    D    D    B
##   X249 X250
## 1    A    C
## 2    A    C
## 3    B    D
## 4    A    B
## 5    D    B
## 6    B    A

Calculating Transition Matrix from Raw Data

F12 = t(dat[1:2,])
F23 = t(dat[2:3,])
F34 = t(dat[3:4,])
F45 = t(dat[4:5,])
F56 = t(dat[5:6,])
#b = t(a)
s = createSequenceMatrix(F34)
s

##    A  B  C  D
## A 12 17 17 14
## B 18 21 25  7
## C  8 19 13 18
## D 16 14 11 20

Calculating Transition Matrix from Raw Data

s = list()
sequenceMatr = list()
for (i in seq(1:5)){
  j = i+1
  s[[i]] = t(dat[c(i:j),])
  sequenceMatr[[i]] <- createSequenceMatrix(s[[i]])
}
sequenceMatr

## [[1]]
##    A  B  C  D
## A 18 15 17 16
## B 19 15 10 15
## C  8 15 19 18
## D  7 15 22 21
## 
## [[2]]
##    A  B  C  D
## A 15 15 16  6
## B 15 17  9 19
## C 11 16 19 22
## D 19 23 14 14
## 
## [[3]]
##    A  B  C  D
## A 12 17 17 14
## B 18 21 25  7
## C  8 19 13 18
## D 16 14 11 20
## 
## [[4]]
##    A  B  C  D
## A  9 13 16 16
## B 21 17 22 11
## C 19  8 20 19
## D 14 13 19 13
## 
## [[5]]
##    A  B  C  D
## A 13 12 24 14
## B 11 18 14  8
## C 17 15 22 23
## D 15 14 13 17

Verifying Time Homogeneity

verifyMarkovProperty(), assessOrder(), assessStationarity() are few other functions in the markov chain package. But these functions expects a sequence data to be passed as input.

verifyHomogeneity(sequenceMatr)

## Testing homogeneity of DTMC underlying input list 
## ChiSq statistic is 74.01893 d.o.f are 60 corresponding p-value is 0.1053601

## $statistic
## [1] 74.01893
## 
## $dof
## [1] 60
## 
## $pvalue
## [1] 0.1053601

Calculating TPM

• The states for CoCo data are:
  • State 1 Purchase Coonoor Lite
  • State 2 Ooty Strong
  • State 3 Both
  • State 4 None

CoCo.data <- as.matrix(read.csv('./CoCo.csv', row.names = 1, header = TRUE))

CoCo.data

##     S1  S2  S3  S4
## S1 520  56  30  29
## S2  30 191  24  21
## S3  38  18 110  17
## S4  11  16  16 123

Markov Chain Object from TPM

CoCo.mc <- new('markovchain', 
                  states = colnames(CoCo.tpm),
                  transitionMatrix = CoCo.tpm,
                  name = 'CoCo MC')

print(CoCo.mc)

##            S1         S2         S3         S4
## S1 0.81889764 0.08818898 0.04724409 0.04566929
## S2 0.11278195 0.71804511 0.09022556 0.07894737
## S3 0.20765027 0.09836066 0.60109290 0.09289617
## S4 0.06626506 0.09638554 0.09638554 0.74096386

Markov Chain Object from TPM

str(CoCo.mc)

## Formal class 'markovchain' [package "markovchain"] with 4 slots
##   ..@ states          : chr [1:4] "S1" "S2" "S3" "S4"
##   ..@ byrow           : logi TRUE
##   ..@ transitionMatrix: num [1:4, 1:4] 0.8189 0.1128 0.2077 0.0663 0.0882 ...
##   .. ..- attr(*, "dimnames")=List of 2
##   .. .. ..$ : chr [1:4] "S1" "S2" "S3" "S4"
##   .. .. ..$ : chr [1:4] "S1" "S2" "S3" "S4"
##   ..@ name            : chr "CoCo MC"

Markov Chain Object from TPM

summary(CoCo.mc)

## CoCo MC  Markov chain that is composed by: 
## Closed classes: 
## S1 S2 S3 S4 
## Recurrent classes: 
## {S1,S2,S3,S4}
## Transient classes: 
## NONE 
## The Markov chain is irreducible 
## The absorbing states are: NONE

Markov Chain Object from TPM

states(CoCo.mc)

## [1] "S1" "S2" "S3" "S4"

CoCo.mc@transitionMatrix

##            S1         S2         S3         S4
## S1 0.81889764 0.08818898 0.04724409 0.04566929
## S2 0.11278195 0.71804511 0.09022556 0.07894737
## S3 0.20765027 0.09836066 0.60109290 0.09289617
## S4 0.06626506 0.09638554 0.09638554 0.74096386

Plotting the Markov Chain Object

plot(CoCo.mc, package = 'diagram' )

Markov Chain Operations: Transition Probability after n Steps

Use the Rain data and CoCo data to calculate steady state probability.

\(P^n\) where P is the TPM

Calculate CoCo TPM after 3 time periods

• Markov Chain Object Method

CoCo.mc ^ 3

## CoCo MC^3 
##  A  4 - dimensional discrete Markov Chain defined by the following states: 
##  S1, S2, S3, S4 
##  The transition matrix  (by rows)  is defined as follows: 
##           S1        S2        S3        S4
## S1 0.6060806 0.1807451 0.1015178 0.1116565
## S2 0.2586340 0.4319802 0.1499508 0.1594350
## S3 0.3634655 0.1922070 0.2763874 0.1679400
## S4 0.1945209 0.1925129 0.1600832 0.4528829

Transition Probability after n Steps / Time Periods

\(P^n\) where P is the TPM Calculate CoCo TPM after 3 time periods - Matrix Exponentiation Method - % ^ % operator from the expm Package

CoCo.tpm %^% 3

##           S1        S2        S3        S4
## S1 0.6060806 0.1807451 0.1015178 0.1116565
## S2 0.2586340 0.4319802 0.1499508 0.1594350
## S3 0.3634655 0.1922070 0.2763874 0.1679400
## S4 0.1945209 0.1925129 0.1600832 0.4528829

Transition Probability after n Steps / Time Periods

\(P^n\) where P is the TPM Calculate CoCo TPM after 3 time periods - Matrix Multiplication Method

CoCo.tpm %*% CoCo.tpm %*% CoCo.tpm 

##           S1        S2        S3        S4
## S1 0.6060806 0.1807451 0.1015178 0.1116565
## S2 0.2586340 0.4319802 0.1499508 0.1594350
## S3 0.3634655 0.1922070 0.2763874 0.1679400
## S4 0.1945209 0.1925129 0.1600832 0.4528829

Transition Probability after n Steps / Time Periods

Exercise

Calculate Rain TPM after 5 time periods:

- Markov Chain Object Method
- TMP Exponentiation Method
- TPM Multiplication Method

Steady State Probability Calculation

Calculate Long Term CoCo Probability

Brute Force Method

for(i in 1:10)  {
    print(paste('i = ', i, sep = ''))
    print(CoCo.tpm %^% i)
    print('')
}

## [1] "i = 1"
##            S1         S2         S3         S4
## S1 0.81889764 0.08818898 0.04724409 0.04566929
## S2 0.11278195 0.71804511 0.09022556 0.07894737
## S3 0.20765027 0.09836066 0.60109290 0.09289617
## S4 0.06626506 0.09638554 0.09638554 0.74096386
## [1] ""
## [1] "i = 2"
##           S1        S2         S3         S4
## S1 0.6933760 0.1445902 0.07944493 0.08258885
## S2 0.1973062 0.5420189 0.13195764 0.12871720
## S3 0.3121105 0.1570176 0.38895141 0.14192049
## S4 0.1342493 0.1559518 0.14118194 0.56861695
## [1] ""
## [1] "i = 3"
##           S1        S2        S3        S4
## S1 0.6060806 0.1807451 0.1015178 0.1116565
## S2 0.2586340 0.4319802 0.1499508 0.1594350
## S3 0.3634655 0.1922070 0.2763874 0.1679400
## S4 0.1945209 0.1925129 0.1600832 0.4528829
## [1] ""
## [1] "i = 4"
##           S1        S2        S3        S4
## S1 0.5451819 0.2039802 0.1167253 0.1341127
## S2 0.3022167 0.3631064 0.1566962 0.1779807
## S3 0.3878390 0.2134396 0.2168351 0.1818863
## S4 0.2442564 0.2147848 0.1664358 0.3745230
## [1] ""
## [1] "i = 5"
##           S1        S2        S3        S4
## S1 0.5025785 0.2189537 0.1272501 0.1512177
## S2 0.3327683 0.3199465 0.1583832 0.1889021
## S3 0.3987512 0.2263216 0.1854500 0.1894771
## S4 0.2836230 0.2282352 0.1670607 0.3210810
## [1] ""
## [1] "i = 6"
##           S1        S2        S3        S4
## S1 0.4726983 0.2286321 0.1345634 0.1641061
## S2 0.3539933 0.2928686 0.1579991 0.1951390
## S3 0.4031259 0.2341785 0.1689942 0.1937015
## S4 0.3139657 0.2362754 0.1653588 0.2844001
## [1] ""
## [1] "i = 7"
##           S1        S2        S3        S4
## S1 0.4516938 0.2349082 0.1396632 0.1737348
## S2 0.3686540 0.2758606 0.1569290 0.1985563
## S3 0.4044572 0.2389944 0.1604254 0.1961230
## S4 0.3369360 0.2410216 0.1629591 0.2590833
## [1] ""
## [1] "i = 8"
##           S1        S2        S3        S4
## S1 0.4368980 0.2389920 0.1432307 0.1808793
## S2 0.3787457 0.2651652 0.1557733 0.2003158
## S3 0.4044718 0.2419604 0.1560056 0.1975622
## S4 0.3541056 0.2437790 0.1605900 0.2415253
## [1] ""
## [1] "i = 9"
##           S1        S2        S3        S4
## S1 0.4264566 0.2416590 0.1457331 0.1861512
## S2 0.3856801 0.2584313 0.1547599 0.2011286
## S3 0.4039959 0.2437954 0.1537559 0.1984528
## S4 0.3668214 0.2453479 0.1585336 0.2292972
## [1] ""
## [1] "i = 10"
##           S1        S2        S3        S4
## S1 0.4190759 0.2434075 0.1474928 0.1900237
## S2 0.3904427 0.2541862 0.1539492 0.2014219
## S3 0.4034049 0.2449356 0.1526326 0.1990269
## S4 0.3761739 0.2462148 0.1568611 0.2207501
## [1] ""

Markovchain Package feature

steadyStates(CoCo.mc)

##             S1        S2        S3      S4
## [1,] 0.4010336 0.2468627 0.1517436 0.20036

Steady State Probability Calculation

Exercise
Calculate Steady State Probability for Rain Markov Chain

- Brute Force Method
- Markovchain package feature

State Probability after n steps given Initial Probability

\(u^{(n)} = uP^n\) where u is Initial Probability and P is TPM

For CoCo MC calculate probability after 4 period if the initial probability was
- S1 40%, S2 25%, S3 15%, S4 20%

initial.prob <- c(0.40, 0.25, 0.15, 0.20)

• Exponentiation Method

initial.prob * CoCo.mc ^ 4

##             S1        S2        S3        S4
## [1,] 0.4006541 0.2473416 0.1516766 0.2003278

State Probability after n steps given Initial Probability

\(u^{(n)} = uP^n\) where u is Initial Probability and P is TPM

• Matrix Multiplication

initial.prob %*% (CoCo.tpm %^% 4) 

##             S1        S2        S3        S4
## [1,] 0.4006541 0.2473416 0.1516766 0.2003278

State Probability after n steps given Initial Probability

Exercise
For Rain Markov Chain calculate probability of 6+ cm of rain after 3 days if the initial probabilities are

- 0cm : 45%  
- 1-5 cm : 35%  
- 6+ cm : 20%  

Absorbing Markov Chains

Airwaves India
- The states for Airwaves India are:

- State 1: Customer Churn that generated no revenue
- State 2: Customer Churn that generated INR 200 profit per month (Incoming calls)
- State 3: Customer State that generated INR 300 profit per month
- State 4: Customer State that generated INR 400 profit per month
- State 5: Customer State that generated INR 600 profit per month
- State 6: Customer State that generated INR 800 profit per month

airwaves.tpm <- as.matrix(read.csv('./Airwaves_India.csv', row.names = 1, header = T))

airwaves.tpm

##      S1   S2  S3   S4   S5   S6
## S1 1.00 0.00 0.0 0.00 0.00 0.00
## S2 0.00 1.00 0.0 0.00 0.00 0.00
## S3 0.05 0.05 0.9 0.00 0.00 0.00
## S4 0.10 0.05 0.0 0.80 0.05 0.00
## S5 0.20 0.10 0.0 0.05 0.60 0.05
## S6 0.10 0.20 0.0 0.00 0.00 0.70

Absorbing Markov Chains

airwaves.mc <- new('markovchain', 
                  states = as.character(c('S1','S2','S3','S4','S5','S6')),
                  transitionMatrix = airwaves.tpm,
                  name = "airwaves's")

Absorbing Markov Chains

summary(airwaves.mc)

## airwaves's  Markov chain that is composed by: 
## Closed classes: 
## S1 
## S2 
## Recurrent classes: 
## {S1},{S2}
## Transient classes: 
## {S3},{S4,S5},{S6}
## The Markov chain is not irreducible 
## The absorbing states are: S1 S2

Absorbing Markov Chains

states(airwaves.mc)

## [1] "S1" "S2" "S3" "S4" "S5" "S6"

transientStates(airwaves.mc)

## [1] "S3" "S4" "S5" "S6"

absorbingStates(airwaves.mc)

## [1] "S1" "S2"

Absorbing Markov Chains

#plotmat(airwaves.mc)
plot(airwaves.mc, package = 'diagram' )

Airwaves’s India - Canonical Form

Arranging the TPM in I, O, R, Q matrix form:

canonicForm(airwaves.mc)

## airwaves's 
##  A  6 - dimensional discrete Markov Chain defined by the following states: 
##  S1, S2, S3, S4, S5, S6 
##  The transition matrix  (by rows)  is defined as follows: 
##      S1   S2  S3   S4   S5   S6
## S1 1.00 0.00 0.0 0.00 0.00 0.00
## S2 0.00 1.00 0.0 0.00 0.00 0.00
## S3 0.05 0.05 0.9 0.00 0.00 0.00
## S4 0.10 0.05 0.0 0.80 0.05 0.00
## S5 0.20 0.10 0.0 0.05 0.60 0.05
## S6 0.10 0.20 0.0 0.00 0.00 0.70

Airwaves’s India - Calculating Q, R & I Matrices

q.matrix <- airwaves.mc[transientStates(airwaves.mc), transientStates(airwaves.mc)]

q.matrix

##     S3   S4   S5   S6
## S3 0.9 0.00 0.00 0.00
## S4 0.0 0.80 0.05 0.00
## S5 0.0 0.05 0.60 0.05
## S6 0.0 0.00 0.00 0.70

r.matrix <- airwaves.mc[transientStates(airwaves.mc), absorbingStates(airwaves.mc)]

r.matrix

##      S1   S2
## S3 0.05 0.05
## S4 0.10 0.05
## S5 0.20 0.10
## S6 0.10 0.20

Absorbing Markov Chains - Fundamental Matrix N

F = (I - Q)^-1
- where I is an Identity Matrix of the same size a Q Matrix

Entry \(\sf{f_{ij}}\) of F
- Expected number of times the process is in transient state \(\sf{s_{j}}\) if it started in transient state \(\sf{s_{i}}\)

f.matrix <- solve(diag(1, nrow = nrow(q.matrix), ncol = ncol(q.matrix)) - q.matrix)

f.matrix

##    S3        S4        S5        S6
## S3 10 0.0000000 0.0000000 0.0000000
## S4  0 5.1612903 0.6451613 0.1075269
## S5  0 0.6451613 2.5806452 0.4301075
## S6  0 0.0000000 0.0000000 3.3333333

Absorbing Markov Chains - Fundamental Matrix N

What is the number of times the airwaves customer is in State 5 given customer started in State 4?

f.matrix['S5', 'S4']

## [1] 0.6451613

Absorbing Markov Chains - Probability of Absorption

B = FR

Entry \(\sf{b_{ij}}\) of B
- Probability an absorbing chain will be absorbed in absobing state \(\sf{s_{j}}\) it started in transient state \(\sf{s_{i}}\)_

absorption.prob <- f.matrix %*% r.matrix

absorption.prob

##           S1        S2
## S3 0.5000000 0.5000000
## S4 0.6559140 0.3440860
## S5 0.6236559 0.3763441
## S6 0.3333333 0.6666667

Absorbing Markov Chains - Probability of Absorption

If the customer is in state 6, calculate the probability of eventual absorption in state 2?

absorption.prob['S6', 'S2']

## [1] 0.6666667

Absorbing Markov Chains - Time to Absorption

t = Fc
- where c is a column vector of all 1s

Entry \(\sf{t_{i}}\) of t
- Expected number of steps before the chain is absorbed if it started in state \(\sf{s_{i}}\)

time.to.absorb <- f.matrix %*% rep(1, nrow(f.matrix))

time.to.absorb

##         [,1]
## S3 10.000000
## S4  5.913978
## S5  3.655914
## S6  3.333333

Absorbing Markov Chains - Time to Absorption

Calculate the expected value of time taken to absorbtion if the current state is S4?

time.to.absorb[2]

## [1] 5.913978

Absorbing Markov Chains - Exercise

ehrenfest.csv give the TPM for Ehrenfest Model. - There are 2 urns with 4 balls between them. - At each step one of the r balls is chosen at random and moved from the urn to the other urn. - States in TPM are based on number of balls in 1st urn. - PS: This model is used to explain diffusion of gases.

1. Make states 0 & 4 absorbing states
2. Calculate
   1. Expected time in different transient states
      
   2. Time to Absorption
      
   3. Absroption Probability