####Question 1#####
# Define the exponential term
exp_term <- exp(-0.7)
# Define matrix P_A
P_A <- matrix(c(
1 - exp_term, exp_term, 0, 0, 0, 0,
1 - exp_term, 0, exp_term, 0, 0, 0,
0, 1 - exp_term, 0, exp_term, 0, 0,
0, 0, 1 - exp_term, 0, exp_term, 0,
0, 0, 0, 1 - exp_term, 0, exp_term,
0, 0, 0, 0, 1 - exp_term, exp_term
), nrow = 6, byrow = TRUE)
# Define matrix P_B
P_B <- matrix(c(
1 - exp_term, exp_term, 0, 0,
1 - exp_term, 0, exp_term, 0,
1 - exp_term, 0, 0, exp_term,
1 - exp_term, 0, 0, exp_term
), nrow = 4, byrow = TRUE)
# Print the matrices to verify
print(P_A)
print(P_B)
####Question 2####
library(expm)
#Take P_A and P_B to the tenth power
P_A_10 <- P_A %^% 10
print(P_A_10)
P_B_10 <- P_B %^% 10
print(P_B_10)
####Question 3####
#Solved algebraically - see Appendix for detailed working
####Question 4#####
# Define the stationary distribution vector
pi <- c(0.17240929, 0.17007036, 0.16776316, 0.16548727, 0.16324224, 0.16102768)
# Define the corresponding discounts for each state
discounts <- c(0, 0.05, 0.10, 0.15, 0.20, 0.25)
# Calculate the average discount
average_discount <- sum(pi * discounts)
# Print the average discount
average_discount
# Define the stationary distribution vector for Scheme (B)
pi_B <- c(0.5034147, 0.24998834, 0.12414054, 0.12245643)
# Define the corresponding discounts for each state in Scheme (B)
discounts_B <- c(0, 0.15, 0.35, 0.55)
# Calculate the average discount for Scheme (B)
average_discount_B <- sum(pi_B * discounts_B)
# Print the average discount for Scheme (B)
average_discount_B
####Question 5####
# Function to simulate multiple antithetic paths
simulate_paths_antithetic_uniform <- function(P, initial_state, n, m) {
paths <- matrix(0, nrow = m, ncol = n + 1)
for (i in 1:(m / 2)) {
# Generate a normal path
path <- numeric(n + 1)
path[1] <- initial_state
# Generate an antithetic path
antithetic_path <- numeric(n + 1)
antithetic_path[1] <- initial_state
for (j in 2:(n + 1)) {
u <- runif(1) # Generate a uniform random variable
cumulative_prob <- cumsum(P[path[j - 1], ]) # Cumulative probabilities
path[j] <- which.max(u <= cumulative_prob)
y <- 1 - u # Generate 1 - u for antithetic path
cumulative_prob_antithetic <- cumsum(P[antithetic_path[j - 1], ]) # Cumulative
probabilities for antithetic path
antithetic_path[j] <- which.max(y <= cumulative_prob_antithetic)
}
paths[i, ] <- path
paths[m - i + 1, ] <- antithetic_path
}
return(paths)
}
# Example usage:
set.seed(1) # Setting seed for reproducibility
initial_state <- 1 # Initial state index
n = 10
m=1000
paths_uni_anti <- simulate_paths_antithetic_uniform(P_A, initial_state, n, m)
# Define the states labels
states_labels <- c(0,5,10,15,20,25)
paths_uni_anti <- apply(paths_uni_anti, 2, function(x) states_labels[x])
paths_uni_anti
end_states_uni_anti <- paths_uni_anti[, n + 1]
proportions_uni_anti <- table(end_states_uni_anti) / m
print(proportions_uni_anti)
df <- as.data.frame(paths_uni_anti)
df$path_id <- 1:m
df_long <- df %>%
pivot_longer(-path_id, names_to = "time", values_to = "state") %>%
mutate(time = as.numeric(gsub("V", "", time)) - 1)
ggplot(df_long, aes(x = time, y = state, group = path_id)) +
geom_line(alpha = 0.1) +
theme_minimal() +
labs(title = "Simulated Markov Chain Paths", x = "Time", y = "State")
####Question 6####
# Simulation for E[X_1 | X_0 = e_1]:
set.seed(12) # Setting seed for reproducibility
initial_state <- 1 # Initial state index
n=1
m=100000
paths_uni_anti_1 <- simulate_paths_antithetic_uniform(P_A, initial_state, n, m)
# Define the states labels
states_labels <- c(0,5,10,15,20,25)
paths_uni_anti_1 <- apply(paths_uni_anti_1, 2, function(x) states_labels[x])
paths_uni_anti_1
end_states_uni_anti_1 <- paths_uni_anti_1[, n + 1]
proportions_uni_anti_1 <- table(end_states_uni_anti_1) / m
print(proportions_uni_anti_1)
# Simulation for E[X_10 | X_0 = e_1]:
set.seed(123) # Setting seed for reproducibility
initial_state <- 1 # Initial state index
n = 10
m=100000
paths_uni_anti_2 <- simulate_paths_antithetic_uniform(P_A, initial_state, n, m)
# Define the states labels
states_labels <- c(0,5,10,15,20,25)
paths_uni_anti_2 <- apply(paths_uni_anti_2, 2, function(x) states_labels[x])
paths_uni_anti_2
end_states_uni_anti_2 <- paths_uni_anti_2[, n + 1]
proportions_uni_anti_2 <- table(end_states_uni_anti_2) / m
print(proportions_uni_anti_2)
(P_A %^% 10)[1,]
####Question 7####
#See technical appendix for detailed working
####Question 8####
# Define the first set of discounts for Scheme A
array1 <- c(0.02897969, 0.04982045, 0.08564885, 0.14724328, 0.25313336,
0.43517437)
array2 <- c(0, 0.05, 0.1, 0.15, 0.2, 0.25)
# Reshape array1 to a column vector (matrix in R with n rows and 1 column)
array1_column <- matrix(array1, ncol = 1)
# Perform matrix multiplication
# t(array1_column) transposes the column vector to a row vector
# %*% is the operator for matrix multiplication in R
result <- t(array1_column) %*% array2
# Calculate leftover by multiplying the result by 3500
leftover <- result * 3500
# Print the results
cat("Scheme A new discount:", result, "\n")
cat("Exp cost of running Scheme A:", leftover, "\n")
# Define the second set of discounts for Scheme B
array1 <- c(0.3677619, 0.23251309, 0.14700363, 0.25272138)
array2 <- c(0, 0.15, 0.35, 0.55)
# Reshape array1 to a column vector (matrix in R with n rows and 1 column)
array1_column <- matrix(array1, ncol = 1)
# Perform matrix multiplication
# t(array1_column) transposes the column vector to a row vector
# %*% is the operator for matrix multiplication in R
result <- t(array1_column) %*% array2
# Calculate leftover by multiplying the result by 3500
leftover <- result * 3500
# Print the results
cat("Scheme B new discount:", result, "\n")
cat("Exp cost of running Scheme B:", leftover, "\n")
