# R commands to find the expected value of a state in a Markov chain
# assuming optimal stopping time.
#
#  AUTHOR: M.V.Wickerhauser
#  DATE: 2018-02-14
#  UPDATE: 2018-02-16
#  TEXT: Gregory F. Lawler, "Introduction to Stochastic Processes,"
#  second edition, ISBN 978-1-58488-651-8, Chapman and Hall/CRC, 2006.

### Preliminary:
# Define a function that returns a vector of maximum coordinates
vecmax <- function(v1,v2) {
 z <- v1;			# initialize output vector z
 ind <- which(v2>v1);		# indices of larger values in v2
 z[ind] <- v2[ind];		# now z = max(v1,v2)
 return(z)
 }
# State space for the dice game:
S <- c(1,2,3,4,5,6);
# Define the dice game transition matrix P:
r1 <- c(1,1,1,1,1,1)/6;  	# constant rows for states 1 through 5
r6 <- c(0,0,0,0,0,1);		# state 6 is absorbing
P <-  matrix(c(r1,r1,r1,r1,r1,r6),nrow=6,byrow=T) # transition matrix
# Define the payoff function: f(k)=k, k=1,...,5, but f(6)=0.
f <- c(1,2,3,4,5,0);
# Set the initial superharmonic function u1:
u1 <- c(5,5,5,5,5,0);


### Example 1, p.91: dice game
# Find the state valuation by iterating:
u <- u1;
for(i in 1:10)  u <- vecmax( f, P %*% u ) # iterate 10 times
u				# publish the result.
for(i in 1:990)  u <- vecmax( f, P %*% u ) # iterate 990 more times
u				# publish the final result.
S2 <- S[u==f]; S2		# Identify the "stop" states  S2
S1 <- S[u>f]; S1		# Identify the "continue" states S1


# Example 1, p.94: dice game with cost
# Define the cost function: cost(k)=1, k=1,...,6
cost=c(1,1,1,1,1,1);
# Find the state valuation by iterating:
u <- u1;
for(i in 1:10)  u <- vecmax( f, P %*% u - cost ) # iterate 10 times
u				# publish the result.
for(i in 1:990)  u <- vecmax( f, P %*% u - cost ) # iterate 990 more times
u				# publish the final result.
S2 <- S[u==f]; S2		# Identify the "stop" states  S2
S1 <- S[u>f]; S1		# Identify the "continue" states S1


# Example 1, p.96: dice game with discounting
# Define the discount factor alpha:
alpha <- 0.8
# Find the state valuation by iterating:
u <- u1;
for(i in 1:10)  u <- vecmax( f, alpha * P %*% u ) # iterate 10 times
u				# publish the result.
for(i in 1:990)  u <- vecmax( f, alpha * P %*% u ) # iterate 990 more times
u				# publish the final result.
S2 <- S[u==f]; S2		# Identify the "stop" states  S2
S1 <- S[u>f]; S1		# Identify the "continue" states S1