Chapter 17: Introduction to Markov Chains

import numpy as np
import matplotlib.pyplot as plt

# Configure plots for better readability
plt.style.use('seaborn-v0_8-whitegrid')
plt.rcParams['figure.figsize'] = (10, 6)

# Define the states
states = ['Free', 'Basic', 'Premium', 'Churned']
state_map = {state: i for i, state in enumerate(states)} # Map state names to indices

# Define the transition matrix P
P = np.array([
    [0.60, 0.20, 0.10, 0.10],  # Transitions from Free
    [0.10, 0.60, 0.20, 0.10],  # Transitions from Basic
    [0.05, 0.10, 0.70, 0.15],  # Transitions from Premium
    [0.00, 0.00, 0.00, 1.00]   # Transitions from Churned (Absorbing)
])

# Verify rows sum to 1
print("Transition Matrix P:")
print(P)
print("\nRow sums:", P.sum(axis=1))

def simulate_path(transition_matrix, state_names, start_state_name, num_steps):
    """Simulates a path through the Markov chain."""
    state_indices = list(range(len(state_names)))
    current_state_index = state_names.index(start_state_name)
    path_indices = [current_state_index]

    for _ in range(num_steps):
        # Get the transition probabilities from the current state
        probabilities = transition_matrix[current_state_index, :]
        
        # Choose the next state based on these probabilities
        next_state_index = np.random.choice(state_indices, p=probabilities)
        path_indices.append(next_state_index)
        
        # Update the current state
        current_state_index = next_state_index
        
        # Optional: Stop if an absorbing state (like Churned) is reached
        if probabilities[current_state_index] == 1.0 and np.sum(probabilities) == 1.0:
           # Check if it's an absorbing state (only loops back to itself)
           if transition_matrix[current_state_index, current_state_index] == 1.0:
              # Fill remaining steps if needed, or break
              # For simplicity here, we just let it stay in the absorbing state
              pass 

    # Convert indices back to state names
    path_names = [state_names[i] for i in path_indices]
    return path_names

# Simulate a path for 12 months starting from 'Free'
start_state = 'Free'
steps = 12
simulated_journey = simulate_path(P, states, start_state, steps)

print(f"Simulated {steps}-month journey starting from {start_state}:")
print(" -> ".join(simulated_journey))

# Calculate the transition matrix after 3 steps (months)
n_steps = 3
P_n = np.linalg.matrix_power(P, n_steps)

print(f"Transition Matrix after {n_steps} steps (P^{n_steps}):")
print(np.round(P_n, 3)) # Round for readability

# Example: Probability of being in 'Premium' after 3 months, starting from 'Free'
start_idx = state_map['Free']
end_idx = state_map['Premium']
prob_free_to_premium_3_steps = P_n[start_idx, end_idx]

print(f"\nProbability(State=Premium at month 3 | State=Free at month 0): {prob_free_to_premium_3_steps:.3f}")

# Example: Probability of being 'Churned' after 6 months, starting from 'Basic'
n_steps_long = 6
P_n_long = np.linalg.matrix_power(P, n_steps_long)
start_idx_basic = state_map['Basic']
end_idx_churned = state_map['Churned']
prob_basic_to_churned_6_steps = P_n_long[start_idx_basic, end_idx_churned]

print(f"Probability(State=Churned at month 6 | State=Basic at month 0): {prob_basic_to_churned_6_steps:.3f}")

# Weather example
W_states = ['Sunny', 'Cloudy', 'Rainy']
W = np.array([
    [0.7, 0.2, 0.1],
    [0.3, 0.5, 0.2],
    [0.2, 0.4, 0.4]
])

# Find eigenvalues and eigenvectors
# We need the *left* eigenvector (v P = lambda v), numpy finds *right* (P u = lambda u).
# The left eigenvector of P for eigenvalue lambda is the right eigenvector of P.T for eigenvalue lambda.
eigenvalues, eigenvectors = np.linalg.eig(W.T)

# Find the eigenvector corresponding to eigenvalue 1
# Due to floating point precision, check for eigenvalues close to 1
idx = np.isclose(eigenvalues, 1)

if not np.any(idx):
    print("No eigenvalue close to 1 found. Check matrix.")
else:
    # Get the eigenvector corresponding to eigenvalue 1
    # Ensure we select only one eigenvector if multiple eigenvalues are close to 1
    # and take the real part as the eigenvector should be real for a stochastic matrix
    stationary_vector_raw = eigenvectors[:, np.where(idx)[0][0]].flatten().real
    
    # Normalize the eigenvector so its components sum to 1
    stationary_distribution = stationary_vector_raw / np.sum(stationary_vector_raw)

    print("Stationary Distribution (Long-run probabilities):")
    for state, prob in zip(W_states, stationary_distribution):
        print(f"  {state}: {prob:.4f}")

    # Verification: pi * W should be equal to pi
    print("\nVerification (pi * W):", np.dot(stationary_distribution, W))

# Simulate many paths to see the distribution of final states
num_simulations = 5000
simulation_length = 24 # Simulate for 2 years
final_states = []

initial_state = 'Basic' # Example starting state

for _ in range(num_simulations):
    path = simulate_path(P, states, initial_state, simulation_length)
    final_states.append(path[-1]) # Get the state after 'simulation_length' months

# Calculate the proportion of simulations ending in each state
from collections import Counter # Efficient way to count
final_state_counts = Counter(final_states)
# Ensure all states are present in the counts, even if 0
for state in states:
    if state not in final_state_counts:
        final_state_counts[state] = 0

final_state_proportions = {state: final_state_counts[state] / num_simulations for state in states}

print(f"Distribution of states after {simulation_length} months starting from '{initial_state}' (based on {num_simulations} simulations):")
for state in states: # Print in consistent order
    print(f"  {state}: {final_state_proportions[state]:.4f}")

# Compare with n-step transition probability calculation
P_long_term = np.linalg.matrix_power(P, simulation_length)
start_idx = state_map[initial_state]
theoretical_probs = P_long_term[start_idx, :]

print(f"\nTheoretical probabilities after {simulation_length} months starting from '{initial_state}':")
for i, state in enumerate(states):
    print(f"  {state}: {theoretical_probs[i]:.4f}")

# Exercise 1 Code/Calculation Space
W_states = ['Sunny', 'Cloudy', 'Rainy']
W = np.array([
    [0.7, 0.2, 0.1],
    [0.3, 0.5, 0.2],
    [0.2, 0.4, 0.4]
])

# Calculate W^2
W_2 = np.linalg.matrix_power(W, 2)

# Probability Sunny -> Rainy in 2 steps
sunny_idx = W_states.index('Sunny')
rainy_idx = W_states.index('Rainy')
prob_sunny_to_rainy_2_steps = W_2[sunny_idx, rainy_idx]
print(f"(Exercise 1) Probability Sunny -> Rainy in 2 steps: {prob_sunny_to_rainy_2_steps:.4f}")

# Exercise 2 Code/Calculation Space
# Gambler's Ruin
gambler_states = [0, 1, 2, 3, 4, 5] # Amount of money
P_gambler = np.array([
    [1.0, 0.0, 0.0, 0.0, 0.0, 0.0], # From 0 (Ruin)
    [0.5, 0.0, 0.5, 0.0, 0.0, 0.0], # From 1
    [0.0, 0.5, 0.0, 0.5, 0.0, 0.0], # From 2
    [0.0, 0.0, 0.5, 0.0, 0.5, 0.0], # From 3
    [0.0, 0.0, 0.0, 0.5, 0.0, 0.5], # From 4
    [0.0, 0.0, 0.0, 0.0, 0.0, 1.0]  # From 5 (Win)
])

def simulate_gambler(transition_matrix, states, start_state_val, max_steps=200): # Increased max_steps slightly
    """Simulates one game until an absorbing state is reached."""
    state_indices = list(range(len(states)))
    current_state_index = states.index(start_state_val)
    #path_indices = [current_state_index] # Path not needed for final state only

    steps = 0
    while steps < max_steps:
        current_state_val = states[current_state_index]
        # Check if in absorbing state
        if current_state_val == 0 or current_state_val == 5:
            return current_state_val # Return final state
            
        probabilities = transition_matrix[current_state_index, :]
        next_state_index = np.random.choice(state_indices, p=probabilities)
        #path_indices.append(next_state_index)
        current_state_index = next_state_index
        steps += 1
        
    # If max_steps reached without hitting absorbing state (unlikely here but good practice)
    return states[current_state_index] 

num_games = 10000 # Increased simulations for better accuracy
start_money = 3
end_results = []
for _ in range(num_games):
    end_results.append(simulate_gambler(P_gambler, gambler_states, start_money))

from collections import Counter
end_counts = Counter(end_results)

ruin_count = end_counts.get(0, 0) # Use .get for safety
win_count = end_counts.get(5, 0)

print(f"\n(Exercise 2) Gambler's Ruin Simulation ({num_games} games starting at ${start_money}):")
print(f"  Proportion ending in Ruin ($0): {ruin_count / num_games:.4f}") # Should be 0.4
print(f"  Proportion ending in Win ($5): {win_count / num_games:.4f}")  # Should be 0.6

# Exercise 3 Code/Calculation Space
# Verification: pi * W = pi

# From previous calculation (ensure these are accurate)
# Recalculate just in case
eigenvalues, eigenvectors = np.linalg.eig(W.T)
idx = np.isclose(eigenvalues, 1)
stationary_vector_raw = eigenvectors[:, np.where(idx)[0][0]].flatten().real
pi_weather = stationary_vector_raw / np.sum(stationary_vector_raw)

result_vector = np.dot(pi_weather, W)

print("\n(Exercise 3) Stationary Distribution Verification:")
print(f"  Calculated pi: {pi_weather}")
print(f"  Result pi * W: {result_vector}")
# Use np.allclose for robust floating point comparison
print(f"  Are they close? {np.allclose(pi_weather, result_vector)}")

# Exercise 4 Code/Calculation Space

# Reload original P and states just in case they were modified
states = ['Free', 'Basic', 'Premium', 'Churned']
state_map = {state: i for i, state in enumerate(states)}
P = np.array([
    [0.60, 0.20, 0.10, 0.10],  # Transitions from Free
    [0.10, 0.60, 0.20, 0.10],  # Transitions from Basic
    [0.05, 0.10, 0.70, 0.15],  # Transitions from Premium
    [0.00, 0.00, 0.00, 1.00]   # Transitions from Churned (Absorbing)
])

P_modified = P.copy() # Start with original subscription matrix

# Modify the 'Churned' row (index 3)
prob_churn_to_free = 0.05
P_modified[3, state_map['Free']] = prob_churn_to_free
P_modified[3, state_map['Churned']] = 1.0 - prob_churn_to_free # Adjust P_33

print("\n(Exercise 4) Modified Transition Matrix:")
print(P_modified)
print("\nIs 'Churned' still absorbing?", P_modified[3, 3] == 1.0) # Should be False
print("The chain is no longer absorbing, as state 'Churned' can transition to 'Free'.")
print("The chain should now be irreducible if all other states can eventually reach Churned and Churned can reach Free.")

# Try finding the stationary distribution for the modified matrix
eigenvalues_mod, eigenvectors_mod = np.linalg.eig(P_modified.T)
idx_mod = np.isclose(eigenvalues_mod, 1)

if not np.any(idx_mod):
    print("\nNo eigenvalue close to 1 found for modified matrix.")
else:
    stationary_vector_raw_mod = eigenvectors_mod[:, np.where(idx_mod)[0][0]].flatten().real
    stationary_distribution_mod = stationary_vector_raw_mod / np.sum(stationary_vector_raw_mod)
    
    print("\nStationary Distribution for Modified Matrix:")
    for state, prob in zip(states, stationary_distribution_mod):
        print(f"  {state}: {prob:.4f}")
    
    print("\nThis represents the long-run proportion of time the system spends in each state.")
    print("Even though customers churn, the small chance of returning means there's a non-zero steady state for all plans.")
    # Verification
    print("\nVerification (pi_mod * P_mod):", np.dot(stationary_distribution_mod, P_modified))