#This code is from https://www.datacamp.com/tutorial/markov-chain-analysis-r
#and me
#Markov Chain Analysis in R

install.packages("markovchain")
install.packages("dplyr")
library(markovchain)
library(dplyr)

#We need to make a markov2 object and a transElec matrix!
#Let's create our own transElec matrix!

transElec <- matrix(c(0.5573215, 0.4426785, 0.0000000,
                      0.0000000, 0.8678118, 0.1321882,
                      0.0000000, 0.0000000, 1.0000000),
                    byrow=TRUE, ncol=3)

print(transElec)

markov2 <- new('markovchain',
               transitionMatrix = transElec, # These are the transition probabilities of a random industry
               states = c('SDR','AE','CW'))

layout <- matrix(c(0,0,0,1,1,0), ncol = 2, byrow = TRUE)
plot(markov2, node.size = 10, layout = layout)

# Extract Q from the transition matrix
Q <- transElec[1:2,1:2]
# Generate It
It <- diag(2)
#Calculate fundamental matrix
N <- solve(It-Q)
#Generate column vector of 1s
one <- t(t(c(1,1)))

# Calculate the expected steps by matrix multiplication
expected <- N%*%one
print(expected)

library(ggplot2)

initState <- c(1,0,0) # The initial state will be SDR Funnel (mimicking setting an initial discovery appointment)

# Initiate probability vectors
SDRProb <- c()
AEProb <- c()
CW <- c()

# Calculate probabilities for 24 steps.
for(k in 1:24){
  nsteps <- initState*markov2^k
  SDRProb[k] <- nsteps[1,1]
  AEProb[k] <- nsteps[1,2]
  CW[k] <- nsteps[1,3]
}

# Make dataframes and merge them
SDRProb <- as.data.frame(SDRProb)
SDRProb$Group <- 'SDR'
SDRProb$Iter <- 1:24
names(SDRProb)[1] <- 'Value'

AEProb <- as.data.frame(AEProb)
AEProb$Group <- 'AE'
AEProb$Iter <- 1:24
names(AEProb)[1] <- 'Value'

CW <- as.data.frame(CW)
CW$Group <- 'CW'
CW$Iter <- 1:24
names(CW)[1] <- 'Value'

steps <- rbind(SDRProb,AEProb,CW)

# Plot the probabilities using ggplot
ggplot(steps, aes(x = Iter, y = Value, col = Group))+
  geom_line() +
  xlab('Chain Step') +
  ylab('Probability') +
  ggtitle('24 Step Chain Probability Prediction')+
  theme(plot.title = element_text(hjust = 0.5))