Markov Chain Model: 123

rm(list=ls())

## data
dat <- c(2,3,3,2,1,1,1,1,1,2,
       3,1,3,2,3,3,2,1,2,2,
       3,2,3,2,3,3,2,2,2,2,
       3,1,3,2,3,3,2,2,1,2,
       3,2,3,2,1,3,2,2,3,2,
       3,1,3,2,3,3,2,2,2,3,
       3,2,3,2,3,3,1,2,3,2,
       3,2,3,2,3,3,1,2,2,2,
       3,2,3,2,1,3,2,1,2,3,
       3,1,3,2,3,3,2,1,2,1)
## Plot data
plot(dat, type = "s")
points(1:100, dat, pch = 19, cex = 0.5)

## Observed transtions
f.m <- matrix(ncol = 3, nrow = 3)
n <- length(dat)
f.m[1,1] <- sum(dat[-n] == 1 & (dat[-1] - dat[-n]) == 0)
f.m[1,2] <- sum(dat[-n] == 1 & (dat[-1] - dat[-n]) == 1)
f.m[1,3] <- sum(dat[-n] == 1 & (dat[-1] - dat[-n]) == 2)
f.m[2,1] <- sum(dat[-n] == 2 & (dat[-1] - dat[-n]) == -1)
f.m[2,2] <- sum(dat[-n] == 2 & (dat[-1] - dat[-n]) == 0)
f.m[2,3] <- sum(dat[-n] == 2 & (dat[-1] - dat[-n]) == 1)
f.m[3,1] <- sum(dat[-n] == 3 & (dat[-1] - dat[-n]) == -2)
f.m[3,2] <- sum(dat[-n] == 3 & (dat[-1] - dat[-n]) == -1)
f.m[3,3] <- sum(dat[-n] == 3 & (dat[-1] - dat[-n]) == 0)
f.m

     [,1] [,2] [,3]
[1,]    4    7    6
[2,]    8   10   24
[3,]    6   24   10

MLE

gamma.m <- matrix(ncol = 3, nrow = 3)
gamma.m[1, ] <- f.m[1, ] / sum(dat[-n] == 1)
gamma.m[2, ] <- f.m[2, ] / sum(dat[-n] == 2)
gamma.m[3, ] <- f.m[3, ] / sum(dat[-n] == 3)

Stationardat Distribution

u.m <- matrix(1, ncol = 3, nrow = 3)
i.m <- diag(3)
delta.m <- c(1, 1, 1) %*% solve(i.m - gamma.m + u.m)
delta.m %*% gamma.m - delta.m  # It is a stationardat distribution

              [,1]         [,2]         [,3]
[1,] -1.665335e-16 5.551115e-17 1.110223e-16

l <- sum(f.m * log(gamma.m))  # Log likelihood

Testing H0: independence

gamma.m_0 <- matrix(ncol = 3, nrow = 3)
gamma.m_0[ ,1] <- sum(dat[-1] == 1) / (length(dat) - 1)
gamma.m_0[ ,2] <- sum(dat[-1] == 2) / (length(dat) - 1)
gamma.m_0[ ,3] <- sum(dat[-1] == 3) / (length(dat) - 1)
gamma.m_0

          [,1]      [,2]      [,3]
[1,] 0.1818182 0.4141414 0.4040404
[2,] 0.1818182 0.4141414 0.4040404
[3,] 0.1818182 0.4141414 0.4040404

l0 <- sum(f.m * log(gamma.m_0))
Q <- 2 * (l - l0)
1 - pchisq(Q, df = 6 - 2)

[1] 0.01366082

So it’s not independece distribution.

test using contengencdat tables

f.m_0 <- gamma.m_0
f.m_0[1, ] <- f.m_0[1, ] * sum(dat[-1] == 1)
f.m_0[2, ] <- f.m_0[2, ] * sum(dat[-1] == 2)
f.m_0[3, ] <- f.m_0[3, ] * sum(dat[-1] == 3)
Q <- sum((f.m - f.m_0)^2 / f.m_0)
1 - pchisq(Q, df = 2 * 2)

[1] 0.01360463

So it’s not independece distribution.