## ----global_options, include=FALSE--------------------------------------------
knitr::opts_chunk$set(fig.width=8.5, fig.height=6, out.width = "70%")
set.seed(123)
library(knitr)
# ANSI colour codes (U+001B) break the LaTeX build, so keep the output plain
options(crayon.enabled = FALSE, cli.num_colors = 1, cli.hyperlink = FALSE)
Sys.setenv(NO_COLOR = "1")
hook_output <- knit_hooks$get("output")
knit_hooks$set(output = function(x, options) {
   x <- gsub("\033\\[[0-9;]*[A-Za-z]", "", x)
   x <- gsub("\033", "", x, fixed = TRUE)
   lines <- options$output.lines
   if (is.null(lines)) {
     return(hook_output(x, options))  # pass to default hook
   }
   x <- unlist(strsplit(x, "\n"))
   more <- "..."
   if (length(lines)==1) {        # first n lines
     if (length(x) > lines) {
       # truncate the output, but add ....
       x <- c(head(x, lines), more)
     }
   } else {
     x <- c(more, x[lines], more)
   }
   # paste these lines together
   x <- paste(c(x, ""), collapse = "\n")
   hook_output(x, options)
 })

## ----load, results='hide', message=FALSE--------------------------------------
library("markovchain")

## ----showClass, echo=FALSE----------------------------------------------------
showClass("markovchain")
showClass("markovchainList")

## ----mcInitLong---------------------------------------------------------------
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))
mcWeather <- new("markovchain", states = weatherStates, byrow = byRow, 
               transitionMatrix = weatherMatrix, name = "Weather")

## ----mcInitShort--------------------------------------------------------------
mcWeather <- new("markovchain", states = c("sunny", "cloudy", "rain"),
                 transitionMatrix = 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), 
                 name = "Weather")

## ----defaultMc----------------------------------------------------------------
defaultMc <- new("markovchain")

## ----intromcList--------------------------------------------------------------
mcList <- new("markovchainList", markovchains = list(mcWeather, defaultMc), 
		          name = "A list of Markov chains")

## ----operations---------------------------------------------------------------
initialState <- c(0, 1, 0)
after2Days <- initialState * (mcWeather * mcWeather)
after7Days <- initialState * (mcWeather ^ 7)
after2Days
round(after7Days, 3)

## ----operations2--------------------------------------------------------------
initialState <- c(0, 1, 0)
after2Days <- (t(mcWeather) * t(mcWeather)) * initialState
after7Days <- (t(mcWeather) ^ 7) * initialState
after2Days
round(after7Days, 3)

## ----fval---------------------------------------------------------------------
fvals<-function(mchain,initialstate,n) {
  out<-data.frame()
  names(initialstate)<-names(mchain)
  for (i in 0:n)
  {
    iteration<-initialstate*mchain^(i)
    out<-rbind(out,iteration)
  }
  out<-cbind(out, i=seq(0,n))
  out<-out[,c(4,1:3)]
  return(out)
}
fvals(mchain=mcWeather,initialstate=c(90,5,5),n=4)

## ----otherMethods-------------------------------------------------------------
states(mcWeather)
names(mcWeather)
dim(mcWeather)

## ----otherMethods2------------------------------------------------------------
name(mcWeather)
name(mcWeather) <- "New Name"
name(mcWeather)

## ----sortMethod---------------------------------------------------------------
markovchain:::sort(mcWeather)

## ----transProb----------------------------------------------------------------
transitionProbability(mcWeather, "cloudy", "rain")
mcWeather[2,3]

## ----printAndShow-------------------------------------------------------------
print(mcWeather)
show(mcWeather)

## ----mcPlot, echo=FALSE, fig.cap="Weather example. Markov chain plot"---------
if (requireNamespace("igraph", quietly = TRUE)) {
  library(igraph)
  plot(mcWeather,layout = layout_with_fr)
  } else {
  message("igraph unavailable")
  }

## ----mcPlotdiagram, echo=FALSE, fig.cap="Weather example. Markov chain plot with diagram"----
if (requireNamespace("diagram", quietly = TRUE)) {
  library(diagram)
  plot(mcWeather, package="diagram", box.size = 0.04)
  } else {
  message("diagram unavailable")
  }

## ----exportImport1------------------------------------------------------------
mcDf <- as(mcWeather, "data.frame")
mcNew <- as(mcDf, "markovchain")
mcDf
mcIgraph <- as(mcWeather, "igraph")

## ----exportImport2------------------------------------------------------------
if (requireNamespace("msm", quietly = TRUE)) {
require(msm)
Q <- rbind ( c(0, 0.25, 0, 0.25),
             c(0.166, 0, 0.166, 0.166),
             c(0, 0.25, 0, 0.25),
             c(0, 0, 0, 0) )
cavmsm <- msm(state ~ years, subject = PTNUM, data = cav, qmatrix = Q, death = 4)
msmMc <- as(cavmsm, "markovchain")
msmMc
  } else {
  message("msm unavailable")
  }

## ----exporImport3-------------------------------------------------------------
if (requireNamespace("etm", quietly = TRUE)) {
library(etm)
data(sir.cont)
sir.cont <- sir.cont[order(sir.cont$id, sir.cont$time), ]
for (i in 2:nrow(sir.cont)) {
  if (sir.cont$id[i]==sir.cont$id[i-1]) {
    if (sir.cont$time[i]==sir.cont$time[i-1]) {
      sir.cont$time[i-1] <- sir.cont$time[i-1] - 0.5
    }
  }
}
tra <- matrix(ncol=3,nrow=3,FALSE)
tra[1, 2:3] <- TRUE
tra[2, c(1, 3)] <- TRUE
tr.prob <- etm::etm(sir.cont, c("0", "1", "2"), tra, "cens", 1)
tr.prob
etm2mc<-as(tr.prob, "markovchain")
etm2mc
  } else {
  message("etm unavailable")
}

## ----fromAndTo, echo=FALSE, fig.cap="The markovchain methods for import and export"----
library(igraph)
importExportGraph<-graph_from_literal(dataframe++markovchain,markovchain-+igraph,
                                 markovchain++matrix,table-+markovchain,msm-+markovchain,etm-+markovchain,
                                 markovchain++sparseMatrix)
plot(importExportGraph,main="Import - Export from and to markovchain objects")

## ----exportImport4------------------------------------------------------------
myMatr<-matrix(c(.1,.8,.1,.2,.6,.2,.3,.4,.3), byrow=TRUE, ncol=3)
myMc<-as(myMatr, "markovchain")
myMc

## ----cchcMcList---------------------------------------------------------------
stateNames = c("H", "I", "D")
Q0 <- new("markovchain", states = stateNames, 
        transitionMatrix =matrix(c(0.7, 0.2, 0.1,0.1, 0.6, 0.3,0, 0, 1), 
        byrow = TRUE, nrow = 3), name = "state t0")
Q1 <- new("markovchain", states = stateNames, 
        transitionMatrix = matrix(c(0.5, 0.3, 0.2,0, 0.4, 0.6,0, 0, 1), 
        byrow = TRUE, nrow = 3), name = "state t1")
Q2 <- new("markovchain", states = stateNames, 
        transitionMatrix = matrix(c(0.3, 0.2, 0.5,0, 0.2, 0.8,0, 0, 1), 
        byrow = TRUE,nrow = 3), name = "state t2")
Q3 <- new("markovchain", states = stateNames, 
          transitionMatrix = matrix(c(0, 0, 1, 0, 0, 1, 0, 0, 1), 
        byrow = TRUE, nrow = 3), name = "state t3")
mcCCRC <- new("markovchainList",markovchains = list(Q0,Q1,Q2,Q3), 
      name = "Continuous Care Health Community")
print(mcCCRC)

## ----cchcMcList2--------------------------------------------------------------
mcCCRC[[1]]
dim(mcCCRC)

## ----conditionalDistr---------------------------------------------------------
conditionalDistribution(mcWeather, "sunny")

## ----steadyStates-------------------------------------------------------------
steadyStates(mcWeather)

## ----gamblerRuin--------------------------------------------------------------
gamblerRuinMarkovChain <- function(moneyMax, prob = 0.5) {
  m <- markovchain:::zeros(moneyMax + 1)
  m[1,1] <- m[moneyMax + 1,moneyMax + 1] <- 1
  states <- as.character(0:moneyMax)
  rownames(m) <- colnames(m) <- states
  
  for(i in 2:moneyMax){ 
    m[i,i-1] <- 1 - prob
    m[i, i + 1] <- prob   
  }
  
  new("markovchain", transitionMatrix = m, 
      name = paste("Gambler ruin", moneyMax, "dim", sep = " "))
}

mcGR4 <- gamblerRuinMarkovChain(moneyMax = 4, prob = 0.5)
steadyStates(mcGR4)

## ----entropy-rate-weather-----------------------------------------------------
entropyRate(mcWeather)
entropyRate(mcWeather, base = exp(1))

## ----topological-entropy------------------------------------------------------
topologicalEntropy(mcWeather)
normalizedEntropyRate(mcWeather)

## ----absorbingStates----------------------------------------------------------
absorbingStates(mcGR4)
absorbingStates(mcWeather)

## ----renaldoMatrix1-----------------------------------------------------------
P <- markovchain:::zeros(10)
P[1, c(1, 3)] <- 1/2;
P[2, 2] <- 1/3; P[2,7] <- 2/3;
P[3, 1] <- 1;
P[4, 5] <- 1;
P[5, c(4, 5, 9)] <- 1/3;
P[6, 6] <- 1;
P[7, 7] <- 1/4; P[7,9] <- 3/4;
P[8, c(3, 4, 8, 10)] <- 1/4;
P[9, 2] <- 1;
P[10, c(2, 5, 10)] <- 1/3;
rownames(P) <- letters[1:10] 
colnames(P) <- letters[1:10]
probMc <- new("markovchain", transitionMatrix = P, 
              name = "Probability MC")
summary(probMc)

## ----summaryDetails-----------------------------------------------------------
summary(mcWeather, details = TRUE)

## ----transientStates----------------------------------------------------------
transientStates(probMc)

## ----probMc2Canonic-----------------------------------------------------------
probMcCanonic <- canonicForm(probMc)
probMc
probMcCanonic

## ----isAccessible-------------------------------------------------------------
is.accessible(object = probMc, from = "a", to = "c")
is.accessible(object = probMc, from = "g", to = "c")

## ----periodicity--------------------------------------------------------------
E <- matrix(0, nrow = 4, ncol = 4)
E[1, 2] <- 1
E[2, 1] <- 1/3; E[2, 3] <- 2/3
E[3,2] <- 1/4; E[3, 4] <- 3/4
E[4, 3] <- 1

mcE <- new("markovchain", states = c("a", "b", "c", "d"), 
		transitionMatrix = E, 
		name = "E")
is.irreducible(mcE)
period(mcE)

## ----mathematica9Mc-----------------------------------------------------------
mathematicaMatr <- markovchain:::zeros(5)
mathematicaMatr[1,] <- c(0, 1/3, 0, 2/3, 0)
mathematicaMatr[2,] <- c(1/2, 0, 0, 0, 1/2)
mathematicaMatr[3,] <- c(0, 0, 1/2, 1/2, 0)
mathematicaMatr[4,] <- c(0, 0, 1/2, 1/2, 0)
mathematicaMatr[5,] <- c(0, 0, 0, 0, 1)
statesNames <- letters[1:5]
mathematicaMc <- new("markovchain", transitionMatrix = mathematicaMatr,
                   name = "Mathematica MC", states = statesNames)

## ----mcMathematics, fig=TRUE, echo=FALSE, fig.align='center', fig.cap="Mathematica 9 example. Markov chain plot."----
plot(mathematicaMc, layout = layout_with_fr)

## ----mathematica9MC, echo=FALSE-----------------------------------------------
summary(mathematicaMc)

## ----fpTime1, eval=FALSE------------------------------------------------------
# .firstpassageKernel <- function(P, i, n){
#   G <- P
#   H <- P[i,]
#   E <- 1 - diag(size(P)[2])
#   for (m in 2:n) {
#     G <- P %*% (G * E)
#     H <- rbind(H, G[i,])
#   }
#   return(H)
# }

## ----fpTime2------------------------------------------------------------------
firstPassagePdF <- firstPassage(object = mcWeather, state = "sunny", 
                                n = 10)
firstPassagePdF[3, 3]

## ----mfpt1--------------------------------------------------------------------
meanFirstPassageTime(mcWeather)

## ----mfpt2--------------------------------------------------------------------
meanFirstPassageTime(mcWeather,"rain")

## ----mfpt3--------------------------------------------------------------------
firstPassagePdF.long <- firstPassage(object = mcWeather, state = "sunny",  n = 100)
sum(firstPassagePdF.long[,"rain"] * 1:100)

## ----kemeny-weather-----------------------------------------------------------
kemenyConstant(mcWeather)

piWeather <- as.numeric(steadyStates(mcWeather))
meanFirstPassageTime(mcWeather) %*% piWeather

## ----sensitivity-weather------------------------------------------------------
S <- sensitivity(mcWeather, "sunny")
S

## ----sensitivity-direction----------------------------------------------------
eps <- 1e-4
predictedChange <- eps * (S["rain", ] - S["sunny", ])
predictedChange

perturbedMatrix <- mcWeather@transitionMatrix
perturbedMatrix["sunny", "sunny"] <- perturbedMatrix["sunny", "sunny"] - eps
perturbedMatrix["sunny", "rain"] <- perturbedMatrix["sunny", "rain"] + eps
mcWeatherPerturbed <- new("markovchain", states = states(mcWeather),
                          transitionMatrix = perturbedMatrix)

as.numeric(steadyStates(mcWeatherPerturbed)) - as.numeric(steadyStates(mcWeather))

## ----spectral-weather---------------------------------------------------------
slem(mcWeather)
spectralGap(mcWeather)

## ----implied-weather----------------------------------------------------------
impliedTimescales(mcWeather)

## ----cycle-walk-odd-----------------------------------------------------------
randomWalkOnCycle <- function(n) {
  cycleP <- matrix(0, n, n)
  for (i in seq_len(n)) {
    cycleP[i, ((i - 2) %% n) + 1] <- 1 / 2
    cycleP[i, (i %% n) + 1] <- 1 / 2
  }
  new("markovchain", transitionMatrix = cycleP,
      states = as.character(seq_len(n)),
      name = paste("Random walk on the", n, "-cycle"))
}

mcCycle7 <- randomWalkOnCycle(7)
slem(mcCycle7)
abs(cos(2 * pi * 3 / 7)) # closed-form SLEM, attained at k = 3 and k = 4

spectralGap(mcCycle7)
impliedTimescales(mcCycle7)

## ----cycle-walk-even----------------------------------------------------------
mcCycle8 <- randomWalkOnCycle(8)
slem(mcCycle8) # 1, attained at k = 4, where cos(2*pi*4/8) = cos(pi) = -1
round(spectralGap(mcCycle8), 12) # numerically 0, up to floating-point error
impliedTimescales(mcCycle8)["tau2"] # the k = 4 mode, sorted first, never decays

## ----reversible-cycle---------------------------------------------------------
is.reversible(mcCycle7)
is.reversible(mcCycle8)

## ----reversible-directed-cycle------------------------------------------------
states3 <- c("a", "b", "c")
directedCycle <- new("markovchain", states = states3,
  transitionMatrix = matrix(c(0, 1, 0,
                              0, 0, 1,
                              1, 0, 0), byrow = TRUE, nrow = 3,
                            dimnames = list(states3, states3)))
is.reversible(directedCycle)

## ----reversible-weather-------------------------------------------------------
is.reversible(mcWeather)

## ----closest-reversible-cycle-------------------------------------------------
approximation <- closestReversible(directedCycle)
approximation$chain
is.reversible(approximation$chain)

## ----closest-reversible-weather-----------------------------------------------
weatherReversible <- closestReversible(mcWeather)
weatherReversible$chain
weatherReversible$distance

steadyStates(mcWeather)
steadyStates(weatherReversible$chain) # unchanged

## ----mixing-weather-----------------------------------------------------------
mixingTime(mcWeather)
mixingTime(mcWeather, epsilon = 0.01) # a tighter threshold needs more steps

## ----mixing-bound-check-------------------------------------------------------
piWeather <- as.numeric(steadyStates(mcWeather))
upperBound <- log(1 / (0.25 * min(piWeather))) / spectralGap(mcWeather)
upperBound
mixingTime(mcWeather) <= upperBound

## ----mixing-periodic-error, error=TRUE----------------------------------------
try({
mixingTime(mcCycle8)
})

## ----lazy-cycle---------------------------------------------------------------
lazyCycle8 <- lazyChain(mcCycle8, alpha = 0.5)
period(mcCycle8)
period(lazyCycle8)
steadyStates(mcCycle8)
steadyStates(lazyCycle8) # unchanged by laziness

## ----lazy-mixing--------------------------------------------------------------
mixingTime(lazyCycle8)

## ----lazy-slower--------------------------------------------------------------
spectralGap(mcCycle7) # aperiodic already, no laziness needed
spectralGap(lazyChain(mcCycle7, alpha = 0.5))
spectralGap(lazyChain(mcCycle7, alpha = 0.9))

## ----redistribute-weather-----------------------------------------------------
redistribute(mcWeather, steps = 5, initial = "rain")
relaxationTime(mcWeather)

## ----timeCorrelations-weather-------------------------------------------------
days <- c("sunny", "sunny", "rain", "cloudy", "sunny")
timeCorrelations(mcWeather, days, timePoints = c(0, 1, 2, 10))
timeRelaxations(mcWeather, days, initial = "rain",
                timePoints = c(0, 1, 2, 10))

## ----autoplot-flow, eval = requireNamespace("ggplot2", quietly = TRUE), fig.cap = "Evolution of the distribution of the weather chain, starting from rain.", fig.height = 3.2----
ggplot2::autoplot(mcWeather, type = "flow", steps = 10, initial = "rain")

## ----autoplot-comparison, eval = requireNamespace("ggplot2", quietly = TRUE), fig.cap = "Transition matrices of the weather chain and of its lazy version with alpha = 0.5.", fig.height = 3.2----
mcLazyWeather <- lazyChain(mcWeather, alpha = 0.5)
mcLazyWeather@name <- "Lazy weather"
ggplot2::autoplot(mcWeather, type = "comparison", other = mcLazyWeather)

## ----mrt-weather--------------------------------------------------------------
meanRecurrenceTime(mcWeather)

## ----mrt-probMc---------------------------------------------------------------
recurrentStates(probMc)
meanRecurrenceTime(probMc)

## ----data-drunkard------------------------------------------------------------
drunkProbs <- markovchain:::zeros(5)
drunkProbs[1,1] <- drunkProbs[5,5] <- 1
drunkProbs[2,1] <- drunkProbs[2,3] <- 1/2
drunkProbs[3,2] <- drunkProbs[3,4] <- 1/2
drunkProbs[4,3] <- drunkProbs[4,5] <- 1/2

drunkMc <- new("markovchain", transitionMatrix = drunkProbs)
drunkMc

## ----rs-drunkard--------------------------------------------------------------
recurrentStates(drunkMc)

## ----ts-drunkard--------------------------------------------------------------
transientStates(drunkMc)

## ----ap-drunkard--------------------------------------------------------------
absorptionProbabilities(drunkMc)

## ----at-drunkard--------------------------------------------------------------
meanAbsorptionTime(drunkMc)

## ----gamblers-ruin-drunkard---------------------------------------------------
ruinMc <- gamblersRuin(upperBound = 4, prob = 0.5, states = as.character(1:5))
identical(unclass(ruinMc@transitionMatrix), unclass(drunkMc@transitionMatrix))

## ----gamblers-ruin-formula----------------------------------------------------
upperBound <- 6
prob <- 0.35
ruinAsymmetric <- gamblersRuin(upperBound = upperBound, prob = prob)
ap <- absorptionProbabilities(ruinAsymmetric)

r <- (1 - prob) / prob
i <- seq_len(upperBound - 1)
theoreticalRuin <- (r^i - r^upperBound) / (1 - r^upperBound)

cbind(state = i, package = as.numeric(ap[as.character(i), "0"]),
      formula = theoreticalRuin)

## ----fundamental-matrix-drunkard----------------------------------------------
fundamentalMatrix(drunkMc)

## ----subchain-submatrix-------------------------------------------------------
Q <- subchain(drunkMc, transientStates(drunkMc), method = "submatrix")
Q
rowSums(Q)

## ----subchain-fundamental-check-----------------------------------------------
solve(diag(nrow(Q)) - Q)

## ----subchain-renormalize-----------------------------------------------------
watchedDrunk <- subchain(drunkMc, transientStates(drunkMc), method = "renormalize")
watchedDrunk
rowSums(watchedDrunk@transitionMatrix)

## -----------------------------------------------------------------------------
committorAB(mcWeather,3,1)

## ----hitting-data-------------------------------------------------------------
M <- markovchain:::zeros(5)
M[1,1] <- M[5,5] <- 1
M[2,1] <- M[2,3] <- 1/2
M[3,2] <- M[3,4] <- 1/2
M[4,2] <- M[4,5] <- 1/2

hittingTest <- new("markovchain", transitionMatrix = M)
hittingProbabilities(hittingTest)

## ----hitting-targets----------------------------------------------------------
hittingProbabilities(hittingTest, targets = c("1", "5"))
all.equal(hittingProbabilities(hittingTest, solver = "direct"),
          hittingProbabilities(hittingTest, solver = "bicgstab"))

## ----hitting-weather----------------------------------------------------------
hittingProbabilities(mcWeather)

## ----stochasticMonotonicity---------------------------------------------------
# Define the state space for income classes
incomeStates <- c("low", "mid", "high")

# Example 1: A stochastically monotone transition matrix
# Higher starting income leads to higher probabilities of reaching/staying in higher brackets
monotoneMatrix <- matrix(c(0.6, 0.3, 0.1,
                           0.3, 0.5, 0.2,
                           0.1, 0.4, 0.5), 
                         byrow = TRUE, nrow = 3, 
                         dimnames = list(incomeStates, incomeStates))

mcIncomeMonotone <- new("markovchain", states = incomeStates, 
                        transitionMatrix = monotoneMatrix, 
                        name = "Monotone Income MC")

is.stochasticallyMonotone(mcIncomeMonotone)

## ----stochasticMonotonicityNonMonotone----------------------------------------
# Example 2: A non-monotone transition matrix
# Here, a person in the "mid" class has an unusually high chance of dropping 
# to "low" (0.7), which is worse than a person already in the "low" class (0.5).
# This violates the stochastic monotonicity property.
notMonotoneMatrix <- matrix(c(0.5, 0.3, 0.2,
                              0.7, 0.2, 0.1,
                              0.1, 0.4, 0.5), 
                            byrow = TRUE, nrow = 3, 
                            dimnames = list(incomeStates, incomeStates))

mcIncomeNotMonotone <- new("markovchain", states = incomeStates, 
                           transitionMatrix = notMonotoneMatrix, 
                           name = "Non-Monotone Income MC")

is.stochasticallyMonotone(mcIncomeNotMonotone)

## ----birth-death--------------------------------------------------------------
bd <- birthDeath(p = c(0.3, 0.4, 0.5), q = c(0.2, 0.3, 0.1))
bd

## ----to-bounded-chain---------------------------------------------------------
reflected <- toBoundedChain(bd, "reflecting")
reflected@transitionMatrix[c(1, 4), ]

## ----urn-model----------------------------------------------------------------
ehrenfest <- urnModel(balls = 6)
ehrenfest@transitionMatrix

piEhrenfest <- as.numeric(steadyStates(ehrenfest))
cbind(package = piEhrenfest, binomial = dbinom(0:6, 6, 0.5))

## ----merge-with---------------------------------------------------------------
drySeason <- new("markovchain", states = c("sunny", "rain"),
                 transitionMatrix = matrix(c(0.9, 0.1, 0.2, 0.8), byrow = TRUE,
                   nrow = 2, dimnames = list(c("sunny", "rain"), c("sunny", "rain"))),
                 name = "Dry season")
wetSeason <- new("markovchain", states = c("sunny", "rain"),
                 transitionMatrix = matrix(c(0.5, 0.5, 0.6, 0.4), byrow = TRUE,
                   nrow = 2, dimnames = list(c("sunny", "rain"), c("sunny", "rain"))),
                 name = "Wet season")

mergeWith(drySeason, wetSeason, gamma = 0.3)

## ----to-nth-order-------------------------------------------------------------
identical(unclass(toNthOrder(mcWeather, 5)@transitionMatrix),
          unclass((mcWeather ^ 5)@transitionMatrix))

## ----popgen-neutral-----------------------------------------------------------
n <- 8
neutralMoran <- populationGeneticsModel(model = "moran", n = n, s = 0)
neutralWF <- populationGeneticsModel(model = "wright-fisher", n = n, s = 0)

i <- 1:(n - 1)
fixation <- function(mc) {
  B <- absorptionProbabilities(mc)
  as.numeric(B[as.character(i), as.character(n)])
}
cbind(i = i, moran = fixation(neutralMoran),
      wrightFisher = fixation(neutralWF), theoretical = i / n)

## ----popgen-selection---------------------------------------------------------
favouredMoran <- populationGeneticsModel(model = "moran", n = n, s = 0.5)
favouredWF <- populationGeneticsModel(model = "wright-fisher", n = n, s = 0.5)

# fixation probability starting from 4 copies out of 8, neutral vs favoured
c(neutral = as.numeric(absorptionProbabilities(neutralMoran)["4", "8"]),
  favoured = as.numeric(absorptionProbabilities(favouredMoran)["4", "8"]))

## ----randomChains-------------------------------------------------------------
fixed <- matrix(NA, 3, 3)
fixed[1, 1] <- 0.5
randomMarkovChain(3, zeros = 2, mask = fixed, seed = 1)@transitionMatrix
dirichletChain(4, diffusion = 1, diagonalBias = 2,
               seed = 1)@transitionMatrix
identityChain(states = c("a", "b"))@transitionMatrix

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

## ----simulatingAListOfMarkovChain---------------------------------------------
patientStates <- rmarkovchain(n = 5, object = mcCCRC, t0 = "H", 
                              include.t0 = TRUE)
patientStates[1:10,]

## ----lumpability-oz-----------------------------------------------------------
states_oz <- c("rainy", "nice", "snowy")

P_oz <- matrix(
  c(0.5,  0.25, 0.25,
    0.5,  0.00, 0.5,
    0.25, 0.25, 0.5),
  byrow = TRUE,
  nrow = 3,
  dimnames = list(states_oz, states_oz)
)

mcOz <- new(
  "markovchain",
  states = states_oz,
  transitionMatrix = P_oz,
  name = "Land of Oz"
)

partition_oz <- list(
  Bad_Weather = c("rainy", "snowy"),
  Nice_Weather = c("nice")
)

is.lumpable(mcOz, partition_oz)

lumped_oz <- lump(mcOz, partition_oz)
lumped_oz

## ----lumpability-nonexact-----------------------------------------------------
P_bad <- matrix(
  c(0.5, 0.3, 0.2,
    0.7, 0.2, 0.1,
    0.1, 0.4, 0.5),
  byrow = TRUE,
  nrow = 3,
  dimnames = list(c("A1", "A2", "B"), c("A1", "A2", "B"))
)

mcBad <- new(
  "markovchain",
  states = c("A1", "A2", "B"),
  transitionMatrix = P_bad
)

partition_bad <- list(
  Group_A = c("A1", "A2"),
  Group_B = c("B")
)

is.lumpable(mcBad, partition_bad)

try(lump(mcBad, partition_bad))

## ----lumpability-forced-------------------------------------------------------
lump(mcBad, partition_bad, force = TRUE)

## ----autolump-example, eval=FALSE---------------------------------------------
# rc <- c("AAA", "AA", "A", "BBB", "BB", "B", "CCC", "D")
# 
# creditMatrix <- matrix(
#   c(90.81, 8.33, 0.68, 0.06, 0.08, 0.02, 0.01, 0.01,
#     0.70, 90.65, 7.79, 0.64, 0.06, 0.13, 0.02, 0.01,
#     0.09, 2.27, 91.05, 5.52, 0.74, 0.26, 0.01, 0.06,
#     0.02, 0.33, 5.95, 85.93, 5.30, 1.17, 1.12, 0.18,
#     0.03, 0.14, 0.67, 7.73, 80.53, 8.84, 1.00, 1.06,
#     0.01, 0.11, 0.24, 0.43, 6.48, 83.46, 4.07, 5.20,
#     0.21, 0.00, 0.22, 1.30, 2.38, 11.24, 64.86, 19.79,
#     0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 100.00),
#   nrow = 8,
#   byrow = TRUE,
#   dimnames = list(rc, rc)
# ) / 100
# 
# creditMc <- new(
#   "markovchain",
#   states = rc,
#   transitionMatrix = creditMatrix,
#   name = "Credit ratings"
# )
# 
# set.seed(123)
# auto_lumped <- autoLump(creditMc, k = 3)
# 
# auto_lumped$partition
# auto_lumped$lumped_chain

## ----aggregate-oz-------------------------------------------------------------
agg_oz <- aggregateStates(mcOz, k = 2)
agg_oz$partition

isTRUE(all.equal(unname(agg_oz$aggregatedChain@transitionMatrix),
                  unname(lumped_oz@transitionMatrix)))

## ----aggregate-oz-kld---------------------------------------------------------
agg_oz$klDivergence

## ----aggregate-dmm------------------------------------------------------------
states_dmm <- as.character(1:4)
P_dmm <- matrix(
  c(0.5, 0.4, 0.0, 0.1,
    0.4, 0.5, 0.1, 0.0,
    0.0, 0.1, 0.5, 0.4,
    0.1, 0.0, 0.4, 0.5),
  byrow = TRUE, nrow = 4,
  dimnames = list(states_dmm, states_dmm)
)
mcDMM <- new("markovchain", states = states_dmm, transitionMatrix = P_dmm,
             name = "Deng-Mehta-Meyn (2009) example")

agg_dmm <- aggregateStates(mcDMM, k = 2)
agg_dmm$partition
agg_dmm$aggregatedChain@transitionMatrix

## ----aggregate-dmm-kld--------------------------------------------------------
agg_dmm$klDivergence

## ----serialize-oz-dictionary--------------------------------------------------
oz_dict <- toDictionary(mcOz)
str(oz_dict)

## ----serialize-oz-roundtrip---------------------------------------------------
identical(fromDictionary(oz_dict)@transitionMatrix, mcOz@transitionMatrix)

## ----serialize-oz-json--------------------------------------------------------
oz_json <- tempfile(fileext = ".json")
toFile(mcOz, oz_json)
cat(readLines(oz_json), sep = "\n")

identical(fromFile(oz_json)@transitionMatrix, mcOz@transitionMatrix)

## ----serialize-oz-csv---------------------------------------------------------
oz_csv <- tempfile(fileext = ".csv")
toFile(mcOz, oz_csv)
cat(readLines(oz_csv), sep = "\n")

fromFile(oz_csv)@name # the package's own default, not "Land of Oz"

## ----serialize-oz-xml---------------------------------------------------------
oz_xml <- tempfile(fileext = ".xml")
toFile(mcOz, oz_xml)
cat(readLines(oz_xml)[1:4], sep = "\n")

## ----export-oz-graph----------------------------------------------------------
cat(toMermaid(mcOz, digits = 2), "\n")

## ----fitMcbyMLE2--------------------------------------------------------------
weatherFittedMLE <- markovchainFit(data = weathersOfDays, method = "mle",name = "Weather MLE")
weatherFittedMLE$estimate
weatherFittedMLE$standardError

## ----fitMcbyLAPLACE-----------------------------------------------------------
weatherFittedLAPLACE <- markovchainFit(data = weathersOfDays, 
                                    method = "laplace", laplacian = 0.01,
                                    name = "Weather LAPLACE")
weatherFittedLAPLACE$estimate

## ----fitSequenceMatrix--------------------------------------------------------
createSequenceMatrix(stringchar = weathersOfDays)

## ----fitSanitizeAbsorbing-----------------------------------------------------
partial <- c("a", "b", "a", "b", "a")
markovchainFit(partial, possibleStates = c("a", "b", "c"),
               sanitize = "uniform")$estimate@transitionMatrix
markovchainFit(partial, possibleStates = c("a", "b", "c"),
               sanitize = "absorbing")$estimate@transitionMatrix

## ----fitSequenceMatrix2-------------------------------------------------------
myMatr<-matrix(c("a","b","b","a","a","b","b","b","b","a","a","a","b","a"),ncol=2)
createSequenceMatrix(stringchar = myMatr,toRowProbs = TRUE)

## ----fitMcbyBootStrap1--------------------------------------------------------
weatherFittedBOOT <- markovchainFit(data = weathersOfDays, 
                                    method = "bootstrap", nboot = 10)
weatherFittedBOOT$estimate
weatherFittedBOOT$standardError

## ----fitMcbyBootStrap2, eval=FALSE--------------------------------------------
# weatherFittedBOOTParallel <- markovchainFit(data = weathersOfDays,
#                                     method = "bootstrap", nboot = 200,
#                                     parallel = TRUE)
# weatherFittedBOOTParallel$estimate
# weatherFittedBOOTParallel$standardError

## ----fitMcbyBootStrap3, eval=FALSE--------------------------------------------
# markovchainFit(data = weathersOfDays, method = "bootstrap", nboot = 200,
#                parallel = TRUE, num.cores = 4)

## ----fitMcbyMLE1--------------------------------------------------------------
weatherFittedMLE$logLikelihood
weatherFittedBOOT$logLikelihood

## ----confint------------------------------------------------------------------
weatherFittedMLE$confidenceInterval
weatherFittedBOOT$confidenceInterval

## ----multinomial--------------------------------------------------------------
multinomialConfidenceIntervals(transitionMatrix = 
        weatherFittedMLE$estimate@transitionMatrix, 
        countsTransitionMatrix = createSequenceMatrix(weathersOfDays))

## ----fitMclists---------------------------------------------------------------
data(holson)
singleMc<-markovchainFit(data=holson[,2:12],name="holson")

## ----fitMclistsFit1, output.lines=20------------------------------------------
mcListFit<-markovchainListFit(data=holson[,2:6],name="holson")
mcListFit$estimate

## ----fitMclistsFit2-----------------------------------------------------------
c1<-c("a","b","a","a","c","c","a")
c2<-c("b")
c3<-c("c","a","a","c")
c4<-c("b","a","b","a","a","c","b")
c5<-c("a","a","c",NA)
c6<-c("b","c","b","c","a")
mylist<-list(c1,c2,c3,c4,c5,c6)
mylistMc<-markovchainFit(data=mylist)
mylistMc

## ----fitAMarkovChainListfromAlist, output.lines=15----------------------------
markovchainListFit(data=mylist)

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

## ----markovchainListPredict---------------------------------------------------
predict(mcCCRC, newdata = c("H", "H"), n.ahead = 5)

## ----markovchainListPredict2--------------------------------------------------
predict(mcCCRC, newdata = c("H", "H"), n.ahead = 5, continue = TRUE)

## ----testIndependence---------------------------------------------------------
data(rain)
assessIndependence(rain$rain)
data(preproglucacon)
assessIndependence(preproglucacon$preproglucacon, method = "G")
chisq.test(createSequenceMatrix(rain$rain), correct = FALSE)$statistic

## ----testIndependencePublished------------------------------------------------
koeberg <- as.character(read.csv(system.file("extdata", "koeberg_wind.csv",
                                             package = "markovchain"))$state)
assessIndependence(koeberg[14:length(koeberg)], method = "G")

## ----test1--------------------------------------------------------------------
sample_sequence <- c("a", "b", "a", "a", "a", "a", "b", "a", "b", "a", 
                     "b", "a", "a", "b", "b", "b", "a")
verifyMarkovProperty(sample_sequence, method = "G")

## ----test1b-------------------------------------------------------------------
burkett_counts <- matrix(
  c(127, 59,
     11,  7,
     36, 20,
      9,  4),
  nrow = 4,
  byrow = TRUE,
  dimnames = list(
    c("previous state 1", "previous state 2",
      "previous state 4", "previous state 5"),
    c("transition to state 1", "transition elsewhere")
  )
)
chisq.test(burkett_counts, correct = FALSE)

## ----test2--------------------------------------------------------------------
data(rain)
assessOrder(rain$rain)

## ----test3--------------------------------------------------------------------
assessStationarity(rain$rain, 10)

## ----test4--------------------------------------------------------------------
structural.zeros <- matrix(FALSE, 3, 3,
                           dimnames = list(c("a", "b", "c"),
                                           c("a", "b", "c")))
structural.zeros["a", "a"] <- TRUE

structural_sequence <- c(
  "a", "b", "a", "c", "a", "b", "c", "a",
  "b", "c", "b", "a", "c", "b", "a", "c",
  "a", "b", "c", "a", "b", "c", "b", "a"
)

assessStationarity(structural_sequence, nblocks = 4,
                   structural.zeros = structural.zeros)

## ----divergence1--------------------------------------------------------------
kullback_counts <- matrix(
  c(51, 11,  8,
    12, 31,  9,
     6, 11, 10),
  nrow = 3,
  byrow = TRUE,
  dimnames = list(c("a", "b", "c"), c("a", "b", "c"))
)

kullback_P <- matrix(
  c(5/8, 1/4, 1/8,
    1/4, 1/2, 1/4,
    1/4, 3/8, 3/8),
  nrow = 3,
  byrow = TRUE,
  dimnames = list(c("a", "b", "c"), c("a", "b", "c"))
)

kullback_mc <- new("markovchain", states = c("a", "b", "c"),
                   transitionMatrix = kullback_P)

verifyEmpiricalToTheoretical(kullback_counts, kullback_mc, method = "G")

## ----divergence2--------------------------------------------------------------
homogeneity_P <- matrix(
  c(.70, .20, .10,
    .20, .50, .30,
    .10, .30, .60),
  nrow = 3,
  byrow = TRUE,
  dimnames = list(c("a", "b", "c"), c("a", "b", "c"))
)

homogeneity_mc <- new("markovchain", states = c("a", "b", "c"),
                      transitionMatrix = homogeneity_P)

set.seed(123)
x1 <- rmarkovchain(500, homogeneity_mc, t0 = "a")
x2 <- rmarkovchain(500, homogeneity_mc, t0 = "a")

verifyHomogeneity(list(x1, x2), method = "G")

## ----pseudobayes--------------------------------------------------------------
pseudoBayesEstimator <- function(raw, apriori){
  v_i <- rowSums(raw) 
  K_i <- numeric(nrow(raw))
  sumSquaredY <- rowSums(raw^2)
  #get numerator
  K_i_num <- v_i^2-sumSquaredY
  #get denominator
  VQ <- matrix(0,nrow= nrow(apriori),ncol=ncol(apriori))
  for (i in 1:nrow(VQ)) {
    VQ[i,]<-v_i[i]*apriori[i,]
  }
  
  K_i_den<-rowSums((raw - VQ)^2)
  
  K_i <- K_i_num/K_i_den
  
  #get the alpha vector
  alpha <- K_i / (v_i+K_i)
  
  #empirical transition matrix
  Emp<-raw/rowSums(raw)
  
  #get the estimate
  out<-matrix(0, nrow= nrow(raw),ncol=ncol(raw))
  for (i in 1:nrow(out)) {
    out[i,]<-alpha[i]*apriori[i,]+(1-alpha[i])*Emp[i,]
  }
  return(out)
}

## ----pseudobayes2-------------------------------------------------------------
trueMc<-as(matrix(c(0.1, .9,.7,.3),nrow = 2, byrow = 2),"markovchain")
aprioriMc<-as(matrix(c(0.5, .5,.5,.5),nrow = 2, byrow = 2),"markovchain")

smallSample<-rmarkovchain(n=20,object = trueMc)
smallSampleRawTransitions<-createSequenceMatrix(stringchar = smallSample)
pseudoBayesEstimator(
  raw = smallSampleRawTransitions, 
  apriori = aprioriMc@transitionMatrix
) - trueMc@transitionMatrix

biggerSample<-rmarkovchain(n=60,object = trueMc)
biggerSampleRawTransitions<-createSequenceMatrix(stringchar = biggerSample)
pseudoBayesEstimator(
  raw = biggerSampleRawTransitions,
  apriori = aprioriMc@transitionMatrix
) - trueMc@transitionMatrix

bigSample<-rmarkovchain(n=100,object = trueMc)
bigSampleRawTransitions<-createSequenceMatrix(stringchar = bigSample)
pseudoBayesEstimator(
  raw = bigSampleRawTransitions,
  apriori = aprioriMc@transitionMatrix
) - trueMc@transitionMatrix

## ----loadAndDoExample---------------------------------------------------------
weatherStates <- c("sunny", "cloudy", "rain")
byRow <- TRUE
weatherMatrix <- matrix(data = c(0.7, 0.2, 0.1, 
                                 0.3, 0.4, 0.3, 
                                 0.2, 0.4, 0.4), 
                        byrow = byRow, nrow = 3, 
                        dimnames = list(weatherStates, weatherStates))
mcWeather <- new("markovchain", states = weatherStates, 
                 byrow = byRow, transitionMatrix = weatherMatrix, 
                 name = "Weather")      
weathersOfDays <- rmarkovchain(n = 365, object = mcWeather, t0 = "sunny")

## ----MAPFit-------------------------------------------------------------------
hyperMatrix<-matrix(c(1, 1, 2, 
                      3, 2, 1,
                      2, 2, 3), 
                    nrow = 3, byrow = TRUE,
                    dimnames = list(weatherStates,weatherStates))
markovchainFit(weathersOfDays[1:200], method = "map", 
               confidencelevel = 0.92, hyperparam = hyperMatrix)
predictiveDistribution(weathersOfDays[1:200], 
                       weathersOfDays[201:365],hyperparam = hyperMatrix) 

## ----MAPFit2------------------------------------------------------------------
hyperMatrix2<- hyperMatrix[c(2,3,1), c(2,3,1)]
markovchainFit(weathersOfDays[1:200], method = "map", 
               confidencelevel = 0.92, hyperparam = hyperMatrix2)
predictiveDistribution(weathersOfDays[1:200], 
                       weathersOfDays[201:365],hyperparam = hyperMatrix2)

## ----inferHyperparam----------------------------------------------------------
inferHyperparam(transMatr = weatherMatrix, scale = c(10, 10, 10))

## ----inferHyperparam2---------------------------------------------------------
inferHyperparam(data = weathersOfDays[1:15])

## ----inferHyperparam3---------------------------------------------------------
hyperMatrix3 <- inferHyperparam(transMatr = weatherMatrix, 
                                scale = c(10, 10, 10))
hyperMatrix3 <- hyperMatrix3$scaledInference
hyperMatrix4 <- inferHyperparam(data = weathersOfDays[1:15])
hyperMatrix4 <- hyperMatrix4$dataInference

## ----MAPandMLE----------------------------------------------------------------
data(preproglucacon)
preproglucacon <- preproglucacon[[2]][1:200]
MLEest <- markovchainFit(preproglucacon, method = "mle")
MAPest <- markovchainFit(preproglucacon, method = "map")
MLEest$estimate
MAPest$estimate

## ----rCtmcInit----------------------------------------------------------------
energyStates <- c("sigma", "sigma_star")
byRow <- TRUE
gen <- matrix(data = c(-3, 3,
                       1, -1), nrow = 2,
              byrow = byRow, dimnames = list(energyStates, energyStates))
molecularCTMC <- new("ctmc", states = energyStates, 
                 byrow = byRow, generator = gen, 
                 name = "Molecular Transition Model")      

## ----rctmcRandom0-------------------------------------------------------------
statesDist <- c(0.8, 0.2)
rctmc(n = 3, ctmc = molecularCTMC, initDist = statesDist, out.type = "df", include.T0 = FALSE)

## ----ctmcRandom1--------------------------------------------------------------
statesDist <- c(0.8, 0.2)
rctmc(n = Inf, ctmc = molecularCTMC, initDist = statesDist, T = 2)

## ----rctmcSteadyStates--------------------------------------------------------
steadyStates(molecularCTMC)

## ----rctmcExpectedTime--------------------------------------------------------
ExpectedTime(molecularCTMC, 1, 2)
probabilityatT(molecularCTMC, 0.5, 1)

## ----rctmcFitting-------------------------------------------------------------
data <- list(c("a", "b", "c", "a", "b", "a", "c", "b", "c"), 
             c(0, 0.8, 2.1, 2.4, 4, 5, 5.9, 8.2, 9))
ctmcFit(data)

## ----mcWeatherQ---------------------------------------------------------------
mcWeatherQ <- expm::logm(mcWeather@transitionMatrix,method='Eigen')
mcWeatherQ

## ----mcWeatherHalfDay---------------------------------------------------------
mcWeatherHalfDayTM <- expm::expm(mcWeatherQ*.5)
mcWeatherHalfDay <- new("markovchain",transitionMatrix=mcWeatherHalfDayTM,name="Half Day Weather Transition Matrix")
mcWeatherHalfDay

## ----ctmcd1-------------------------------------------------------------------
if(requireNamespace(package='ctmcd', quietly = TRUE)) {
require(ctmcd)
require(expm)
#defines a function to transform a GM into a TM
gm_to_markovchain<-function(object, t=1) {
  if(!inherits(object, c("gm","matrix","Matrix")))
    stop("Error! Expecting either a matrix or a gm object")
  if (inherits(object, c("matrix","Matrix"))) generator_matrix<-object else generator_matrix<-as.matrix(object[["par"]])
  #must add importClassesFrom("markovchain",markovchain) in the NAMESPACE
  #must add importFrom(expm, "expm")
  transitionMatrix<-expm(generator_matrix*t)
  out<-as(transitionMatrix,"markovchain")
  return(out)
}
#loading ctmcd dataset
data(tm_abs)
gm0=matrix(1,8,8) #initializing
diag(gm0)=0
diag(gm0)=-rowSums(gm0)
gm0[8,]=0
gmem=gm(tm_abs,te=1,method="EM",gmguess=gm0) #estimating GM
mc_at_2=gm_to_markovchain(object=gmem, t=2) #converting to TM at time 2
} else {
  warning('package ctmcd unavailable')
}

## ----weatPred1----------------------------------------------------------------

mcWP <- new("markovchain", states = c("rainy", "nice", "snowy"),
         transitionMatrix = matrix(c(0.5, 0.25, 0.25,
                                   0.5, 0, 0.5,
                                   0.25,0.25,0.5), byrow = T, nrow = 3))

## ----weatPred2----------------------------------------------------------------
W0 <- t(as.matrix(c(0, 1, 0)))
W1 <- W0 * mcWP; W1

W2 <- W0 * (mcWP ^ 2); W2

W3 <- W0 * (mcWP ^ 3); W3

## ----weatPred3----------------------------------------------------------------
W7 <- W0 * (mcWP ^ 7)
W7

## ----weatPred4----------------------------------------------------------------
q <- steadyStates(mcWP)
q

## ----weatPred5----------------------------------------------------------------
R0 <- t(as.matrix(c(1, 0, 0)))
R7 <- R0 * (mcWP ^ 7); R7

S0 <- t(as.matrix(c(0, 0, 1)))
S7 <- S0 * (mcWP ^ 7); S7

## ----Alofi1-------------------------------------------------------------------
data("rain", package = "markovchain")
table(rain$rain)

## ----Alofi2-------------------------------------------------------------------
mcAlofi <- markovchainFit(data = rain$rain, name = "Alofi MC")$estimate
mcAlofi

## ----Alofi3-------------------------------------------------------------------
steadyStates(mcAlofi)

## ----ratings1-----------------------------------------------------------------
rc <- c("AAA", "AA", "A", "BBB", "BB", "B", "CCC", "D")
creditMatrix <- matrix(
  c(90.81, 8.33, 0.68, 0.06, 0.08, 0.02, 0.01, 0.01,
    0.70, 90.65, 7.79, 0.64, 0.06, 0.13, 0.02, 0.01,
    0.09, 2.27, 91.05, 5.52, 0.74, 0.26, 0.01, 0.06,
    0.02, 0.33, 5.95, 85.93, 5.30, 1.17, 1.12, 0.18,
    0.03, 0.14, 0.67, 7.73, 80.53, 8.84, 1.00, 1.06,
    0.01, 0.11, 0.24, 0.43, 6.48, 83.46, 4.07, 5.20,
    0.21, 0, 0.22, 1.30, 2.38, 11.24, 64.86, 19.79,
    0, 0, 0, 0, 0, 0, 0, 100
   )/100, 8, 8, dimnames = list(rc, rc), byrow = TRUE)

## ----ratings2-----------------------------------------------------------------
creditMc <- new("markovchain", transitionMatrix = creditMatrix, 
                name = "S&P Matrix")
absorbingStates(creditMc)

## ----economicAnalysis1--------------------------------------------------------
statesNames <- c("customer", "non customer")
P <- markovchain:::zeros(2); P[1, 1] <- .9; P[1, 2] <- .1; P[2, 2] <- .95; P[2, 1] <- .05;
rownames(P) <- statesNames; colnames(P) <- statesNames
mcP <- new("markovchain", transitionMatrix = P, name = "Telephone company")
M <- markovchain:::zeros(2); M[1, 1] <- -20; M[1, 2] <- -30; M[2, 1] <- -40; M[2, 2] <- 0

## ----economicAnalysis2--------------------------------------------------------
c1 <- 100 + conditionalDistribution(mcP, state = "customer") %*% M[1,]
c2 <- 0 + conditionalDistribution(mcP, state = "non customer") %*% M[2,]

## ----economicAnalysis3--------------------------------------------------------
as.numeric((c(1, 0)* mcP ^ 5) %*% (as.vector(c(c1, c2))))

## ----tauchen-rouwenhorst------------------------------------------------------
rho <- 0.979
sigmaEps <- 0.0072
theoreticalVar <- sigmaEps^2 / (1 - rho^2)
theoreticalVar

relativeError <- function(size) {
  outT <- tauchen(alpha = 0, sigma = sigmaEps, rho = rho, size = size)
  outR <- rouwenhorst(alpha = 0, sigma = sigmaEps, rho = rho, size = size)
  piT <- as.numeric(steadyStates(outT$chain))
  piR <- as.numeric(steadyStates(outR$chain))
  varT <- sum(piT * (outT$states - sum(piT * outT$states)) ^ 2)
  varR <- sum(piR * (outR$states - sum(piR * outR$states)) ^ 2)
  c(tauchen = 100 * abs(varT - theoreticalVar) / theoreticalVar,
    rouwenhorst = 100 * abs(varR - theoreticalVar) / theoreticalVar)
}

sapply(c(5, 10, 25), relativeError)

## ----bonusMalus1--------------------------------------------------------------

getBonusMalusMarkovChain <- function(lambda) {
	bmMatr <- markovchain:::zeros(5)
	bmMatr[1, 1] <- dpois(x = 0, lambda)
	bmMatr[1, 3] <- dpois(x = 1, lambda)
	bmMatr[1, 5] <- 1 - ppois(q = 1, lambda)
	
	bmMatr[2, 1] <- dpois(x = 0, lambda)
	bmMatr[2, 4] <- dpois(x = 1, lambda)
	bmMatr[2, 5] <- 1 - ppois(q = 1, lambda)
	
	bmMatr[3, 2] <- dpois(x = 0, lambda)
	bmMatr[3, 5] <- 1 - dpois(x=0, lambda)
 
	bmMatr[4, 3] <- dpois(x = 0, lambda)
	bmMatr[4, 5] <- 1 - dpois(x = 0, lambda)
  
	bmMatr[5, 4] <- dpois(x = 0, lambda)
	bmMatr[5, 5] <- 1 - dpois(x = 0, lambda)
	stateNames <- as.character(1:5)
	out <- new("markovchain", transitionMatrix = bmMatr, 
             states = stateNames, name = "BM Matrix")
	return(out)
}

## ----bonusMalus2--------------------------------------------------------------
bmMc <- getBonusMalusMarkovChain(0.05)
as.numeric(steadyStates(bmMc))

## ----bonusMalus3--------------------------------------------------------------
sum(as.numeric(steadyStates(bmMc)) * c(0.5, 0.7, 0.9, 1, 1.25))

## ----healthIns6---------------------------------------------------------------
ltcDemoPath<-system.file("extdata", "ltdItaData.txt", 
                         package = "markovchain")
ltcDemo<-read.table(file = ltcDemoPath, header=TRUE, 
                    sep = ";", dec = ".")
head(ltcDemo)

## ----healthIns7---------------------------------------------------------------
ltcDemo<-transform(ltcDemo,
                   pIA=0,
                   pII=1-pID,
                   pDD=1,
                   pDA=0,
                   pDI=0)

## ----healthIns8---------------------------------------------------------------
possibleStates<-c("A","I","D")
getMc4Age<-function(age) {
  transitionsAtAge<-ltcDemo[ltcDemo$age==age,]
  
  myTransMatr<-matrix(0, nrow=3,ncol = 3,
                      dimnames = list(possibleStates, possibleStates))
  myTransMatr[1,1]<-transitionsAtAge$pAA[1]
  myTransMatr[1,2]<-transitionsAtAge$pAI[1]
  myTransMatr[1,3]<-transitionsAtAge$pAD[1]
  myTransMatr[2,2]<-transitionsAtAge$pII[1]
  myTransMatr[2,3]<-transitionsAtAge$pID[1]
  myTransMatr[3,3]<-1
  
  myMc<-new("markovchain", transitionMatrix = myTransMatr,
            states = possibleStates,
            name = paste("Age",age,"transition matrix"))
  
  return(myMc)

}

## ----healthIns8-prob----------------------------------------------------------
getFullTransitionTable<-function(age){
  ageSequence<-seq(from=age, to=120)
  k=1
  myList=list()
  for ( i in ageSequence) {
    mc_age_i<-getMc4Age(age = i)
    myList[[k]]<-mc_age_i
    k=k+1
  }
  myMarkovChainList<-new("markovchainList", markovchains = myList,
                         name = paste("TransitionsSinceAge", age, sep = ""))
  return(myMarkovChainList)
}
transitionsSince100<-getFullTransitionTable(age=100)

## ----healthIns9---------------------------------------------------------------
rmarkovchain(n = 10, object = transitionsSince100,
             what = "matrix", t0 = "A", include.t0 = TRUE)

## ----healthIns10--------------------------------------------------------------
transitionsSince80<-getFullTransitionTable(age=80)
lifeTrajectories<-rmarkovchain(n=1e2, object=transitionsSince80,
                               what="matrix",t0="A",include.t0=TRUE)
temp<-matrix(0,nrow=nrow(lifeTrajectories),ncol = ncol(lifeTrajectories))
temp[lifeTrajectories=="I"]<-1
expected_period_disabled<-mean(rowSums((temp)))
expected_period_disabled

## ----healthIns11--------------------------------------------------------------
mean(rowMeans(12000*temp%*%( matrix((1+0.02)^-seq(from=0, to=ncol(temp)-1)))))

## ----blandenEtAlii------------------------------------------------------------
data("blanden")
mobilityMc <- as(blanden, "markovchain")
mobilityMc

## ----mobility, fig=TRUE, echo=FALSE, fig.align='center', fig.cap="1970 UK cohort mobility data."----
plot(mobilityMc, main = '1970 mobility',vertex.label.cex = 2,
		layout = layout_with_fr)

## ----blandenEtAlii3-----------------------------------------------------------
round(steadyStates(mobilityMc), 2)

## ----preproglucacon1----------------------------------------------------------
data("preproglucacon", package = "markovchain")

## ----preproglucacon2----------------------------------------------------------
mcProtein <- markovchainFit(preproglucacon$preproglucacon, 
                          name = "Preproglucacon MC")$estimate
mcProtein

## ----epid1--------------------------------------------------------------------
craigSendiMatr <- matrix(c(682, 33, 25,
              154, 64, 47,
              19, 19, 43), byrow = T, nrow = 3)
hivStates <- c("0-49", "50-74", "75-UP")
rownames(craigSendiMatr) <- hivStates
colnames(craigSendiMatr) <- hivStates
craigSendiTable <- as.table(craigSendiMatr)
mcM6 <- as(craigSendiTable, "markovchain")
mcM6@name <- "Zero-Six month CD4 cells transition"
mcM6

## ----epid2--------------------------------------------------------------------
eig <- eigen(mcM6@transitionMatrix)
D <- diag(eig$values)

## ----epid3--------------------------------------------------------------------
V <- eig$vectors 
V %*% D %*% solve(V)
d <- D ^ (1/6)
M <- V %*% d %*% solve(V)
mcM1 <- new("markovchain", transitionMatrix = M, states = hivStates)