A markoov proecss is depends only on the state in which it currently is and not on the previous or future states. A markov chain is a particular kind of markov process: discrete time.

A markov chain is defined by some States belonging to a finite State Space S, and by a Transition matrix T. The transition matrix defines the probability to be in a certain state j tomorrow, known that you are in a state i today.

Regular markov chain

A markov chain is REGULAR, if the transition matrix T is regular, so if for some integer \(n\) all entries of \(T^n\) are positive.

Regularity ensures that the MC has a unique stationary distribution and that all states communicate with each other.

##      [,1] [,2]
## [1,]  0.3 0.52
## [2,]  0.7 0.48

This is regular because the entries are all positive, in this case it is abvious because I already have all positive values.

T = matrix(c(0.4,0.6, 1,0),
           nrow = 2,
           ncol = 2)
print(T)

##      [,1] [,2]
## [1,]  0.4    1
## [2,]  0.6    0

In this example, we don’t know if the matrix is regular or not.

To say that it is regular I have to find at least a power of the matrix for which the entries are all positive.

library(expm)

## Caricamento del pacchetto richiesto: Matrix

## 
## Caricamento pacchetto: 'expm'

## Il seguente oggetto è mascherato da 'package:Matrix':
## 
##     expm

m2 = T %^% 2

print(m2)

##      [,1] [,2]
## [1,] 0.76  0.4
## [2,] 0.24  0.6

In this case we have found a positive matrix, so the markov chain is regular.

Another example:

T = matrix(c(1,0, 0.5,0.5),
           nrow = 2,
           ncol = 2)
paste('T:')

## [1] "T:"

print(T)

##      [,1] [,2]
## [1,]    1  0.5
## [2,]    0  0.5

m2 = T %^% 2
paste('T^2:')

## [1] "T^2:"

print(m2)

##      [,1] [,2]
## [1,]    1 0.75
## [2,]    0 0.25

m3 = T %^% 3
paste('T^3:')

## [1] "T^3:"

print(m3)

##      [,1]  [,2]
## [1,]    1 0.875
## [2,]    0 0.125

This markov chain is not regular!

Long run behaviour of MC (steady state)

If the MC is regular, I have that the long run markov chain state \(\pi\) is positive and does not depends on the initial state.

How to comput the long run state ?

I should multiply the Transition matrix by itself “infinite” times !

Or it could be demonstrated that is the same as computing the left eigenvectors of the transition matrix corresponding to eginevalues of 1.

Mathematically, is the solution of the following equation:

\[\pi = \pi \cdot P\]

Interpretation:

\(\pi_{j}\) is the long run probability of finding the process {Xn} in state j.
\(\pi_{j}\) is the long run mean fraction of time that the chain spans in state j.

# Install and load the markovchain package
# install.packages("markovchain")
library(markovchain)

# Step 2: Create the transition matrix
transition_matrix <- matrix(c(0.7, 0.3, 0.2, 0.8), nrow = 2, byrow = TRUE)
print(transition_matrix)

##      [,1] [,2]
## [1,]  0.7  0.3
## [2,]  0.2  0.8

rownames(transition_matrix) <- colnames(transition_matrix) <- c("A", "B")

# Step 3: Create a Markov chain object
mc <- new("markovchain", states = c("A", "B"), transitionMatrix = transition_matrix)

# Step 4: Compute the steady state
steady_state <- steadyStates(mc)

# Print the steady state vector
print(steady_state)

##        A   B
## [1,] 0.4 0.6

It is the same of multipling the same transition matrix infinite times (in this case 2000 times):

mc_2 = transition_matrix %^% 200

print(mc_2)

##     A   B
## A 0.4 0.6
## B 0.4 0.6

Reachability, Communicate, Reducibility and Periodicity

1. Reachability

If exists a path that connect the state i and j. The path can be indirect, but it should exists.

2. Communicate ( i ⟷ j )

If i is reachable from j and if j is reachable from i

Note: Equivalence and Classes

” ⟷ “: is an equivalence relation:

reflexive: i ⟷ i pii = 1
symmetric: i ⟷ j <=> j ⟷ i
transitive: if i ⟷ k & k ⟷ j => i ⟷ j

Hence:

All the states within a Class communicate, so If I am in a state I can come to every other state inside a class.

The state space of the MC is a union of classes.

S = C0 ∪ C1 ∪ C2 ∪… … ∪ Cn

where:

Ci is non empty
Ci ⋂ Cj is empty
Ci are equivalent classes induced by “⟷”

So I can split the MC in different classes.

3. Irreducibility

If there is only one communication class in the markov chain, it is irruducible.

N.B.: if Regulare => Irreducible

4. Period

For an irriducible MC, the PERIOD is the greatest common divisor of n such that Pii(n) >0.

If the greatest common divisor is 1 => APERIODIC.

N.B.: All the states in the same class have the same period.

Recurrence and Transience

1. Recurrence (Loop)

If the probability of returning to i given that you started in i is 1.

3. Transient

If it is not recurrent

Note: Recurrence and Transience are properties of the Class, so all elements inside a class have the same propertiy.

Further details

If \(\mu_{i}\) is the mean time it takes to come back to a certain state i, I can define:

if i is transient => \(\mu\) = ∞ (I’m not coming back)
if i is recurrent:
1. \(\mu\) = ∞ : NULL RECURRENT (I will come back in an infinite amount of time)
2. \(\mu\) < ∞ : POSITIVE RECURRENT (I will come back in a finite amount of time)

Theorems

1. Stationary distribution

if MC is irreducible and positive recurrent:

\[ \pi_{i}=\frac{1}{\mu_{i}} \]

\(\pi_{i}\) is called STATIONARY DISTRIBUTION of the chain.

In fact if I start in i with probability \(\pi_{i}\), then the next day I’ll be still in i with \(\pi_{i}\) probability, and so on…

For this reason is called stationary, it doesn’t change.

2. Theorem

if MC is irreducible and aperiodic:

\[ \pi_{i}=\lim_{n \to \infty} p_{i,i}^{n} = \frac{1}{\mu_{i}} \]

if MC is irreducible and periodic:

\[ \pi_{i}=\lim_{n \to \infty} p_{i,i}^{n\cdot Period} = \frac{1}{\mu_{i}} \]

Markov Chain with ‘markovchain’ package

states: a character vector, listing the states for which transition probabilities are defined.
byrow: a logical element, indicating whether transition probabilities are shown by row or by column.
transitionMatrix: the probabilities of the transition matrix.
name: optional character element to name the DTMC.

library('markovchain')

weatherStates <- c("sunny", "cloudy", "rain")
byRow <- TRUE

weatherMatrix <- matrix(data = c(0.70, 0.2, 0.1, 
                                 0.3, 0.4, 0.3, 
                                 0.2, 0.45, 0.35), 
                        byrow = byRow, 
                        nrow = 3, 
                        dimnames = list(weatherStates, weatherStates))

To initialize the Markov chain I use new()

mcWeather <- new("markovchain", states = weatherStates, byrow = byRow,   
                 transitionMatrix = weatherMatrix, name = "Weather")
print(mcWeather)

##        sunny cloudy rain
## sunny    0.7   0.20 0.10
## cloudy   0.3   0.40 0.30
## rain     0.2   0.45 0.35

What is the probability distribution of expected weather states in two and seven days, given the actual state to be cloudy?

initialState <- c(0, 1, 0)

after2Days <- initialState * (mcWeather * mcWeather)
after7Days <- initialState * (mcWeather ^ 7)

after2Days

##      sunny cloudy  rain
## [1,]  0.39  0.355 0.255

To access the single transition probability:

transitionProbability(mcWeather, "cloudy", "rain")

## [1] 0.3

A function I personally like is that I can also plot the markov chain with plot:

plot(mcWeather)

The plot function also uses communicatingClasses function to separate out states of different communicating classes. All states that belong to one class have same color.

Conditional distributions

conditionalDistribution(mcWeather, "sunny")

##  sunny cloudy   rain 
##    0.7    0.2    0.1

Classification of states

To know all the states of the MC I can use the summary method:

summary(mcWeather)

## Weather  Markov chain that is composed by: 
## Closed classes: 
## sunny cloudy rain 
## Recurrent classes: 
## {sunny,cloudy,rain}
## Transient classes: 
## NONE 
## The Markov chain is irreducible 
## The absorbing states are: NONE

We can easily get the steady states of the MC:

steadyStates(mcWeather)  # 'PI GRECO' IN THEORY

##          sunny    cloudy      rain
## [1,] 0.4636364 0.3181818 0.2181818

The function is.accessible permits to investigate whether a state sj is accessible from state si, that is whether the probability to eventually reach sj starting from si is greater than zero.

is.accessible(object = mcWeather, from = "sunny", to = "rain")

## [1] TRUE

The same could be done to understand if it is irreducible:

is.irreducible(mcWeather)

## [1] TRUE

If the MC is irreducible all the states share the same periodicity.

If the MC is irreducible I can return the period:

period(mcWeather)

## [1] 1

Simulation

weathersOfDays <- rmarkovchain(n = 365, object = mcWeather, t0 = "sunny")
weathersOfDays[1:30]

##  [1] "sunny"  "sunny"  "sunny"  "cloudy" "sunny"  "sunny"  "sunny"  "sunny" 
##  [9] "cloudy" "sunny"  "sunny"  "sunny"  "cloudy" "sunny"  "rain"   "rain"  
## [17] "cloudy" "sunny"  "cloudy" "cloudy" "cloudy" "cloudy" "sunny"  "sunny" 
## [25] "sunny"  "sunny"  "rain"   "cloudy" "rain"   "cloudy"

Estimation

markovchainFit(data, method) will estimate the transition matrix that best fits to your data.

weatherFittedMLE <- markovchainFit(data = weathersOfDays, method = "mle")
print(weatherFittedMLE$estimate)

##           cloudy      rain     sunny
## cloudy 0.3589744 0.3076923 0.3333333
## rain   0.4719101 0.3820225 0.1460674
## sunny  0.2151899 0.1202532 0.6645570

print(weatherFittedMLE$standardError)

##            cloudy       rain      sunny
## cloudy 0.05539095 0.05128205 0.05337605
## rain   0.07281731 0.06551631 0.04051181
## sunny  0.03690476 0.02758797 0.06485412

Prediction

The 3-days forward predictions:

predict(object = weatherFittedMLE$estimate, 
        newdata = c("cloudy", "sunny"),
        n.ahead = 3)

## [1] "sunny" "sunny" "sunny"

Markov Chains