Why Bellman Allocation Problems Need Backward Induction

# Python implementation of bellman equation in Merton Model


# Python implementation of bellman equation in Merton Model












Simple code to find optimal allocation in year t between investing in stocks (X) vs prepaying mortgage (1-X) from simulated joint distribution of stocks and mortgage rates. Can't get the simulation on the Bellman equation to work...

```
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

# Constants
START_MORTGAGE = 500000
INITIAL_EARNINGS = 220000
EARNINGS_GROWTH = 0.05
NON_MORTGAGE_SPENDING = 70000
SPENDING_GROWTH = 0.06
INITIAL_STOCKS = 60000
PROPERTY_VALUE = 650000
PROPERTY_GROWTH = 0.02
RISK_AVERSION = 2
DISCOUNT_RATE = 0.03
MORTGAGE_TERM = 25
SIMULATIONS = 1000
YEARS = 10

# Distribution parameters
STOCK_MEAN_RETURN = 0.09
STOCK_STD_DEV = 0.20
MORTGAGE_MEAN_RATE = 0.03
MORTGAGE_STD_DEV = 0.01
CORRELATION = 0.2

# Generate joint distribution
np.random.seed(42)
covariance = STOCK_STD_DEV * MORTGAGE_STD_DEV * CORRELATION
mean = [STOCK_MEAN_RETURN, MORTGAGE_MEAN_RATE]
cov = [[STOCK_STD_DEV**2, covariance], [covariance, MORTGAGE_STD_DEV**2]]
simulated_data = np.random.multivariate_normal(mean, cov, (SIMULATIONS, YEARS))
stock_returns = simulated_data[:, :, 0]
mortgage_rates = simulated_data[:, :, 1]

def utility(wealth, gamma):
    """CRRA Utility Function."""
    if wealth <= 0:
        return -np.inf
    return (wealth ** (1 - gamma)) / (1 - gamma)

def calculate_excess_earnings(gross_income, tax, spending, mortgage_payment):
    """Calculate excess earnings."""
    return gross_income - tax - spending - mortgage_payment

def tax_calculation(gross_income):
    """Simple progressive tax calculation."""
    personal_allowance = 12570
    basic_rate = 50270
    higher_rate = 150000

    if gross_income <= personal_allowance:
        return 0
    elif gross_income <= basic_rate:
        return (gross_income - personal_allowance) * 0.2
    elif gross_income <= higher_rate:
        return (basic_rate - personal_allowance) * 0.2 + (gross_income - basic_rate) * 0.4
    else:
        return (basic_rate - personal_allowance) * 0.2 + (higher_rate - basic_rate) * 0.4 + (gross_income - higher_rate) * 0.45

def calculate_pmt(rate, term, principal):
    """Calculate the mortgage payment using the PMT formula."""
    if rate == 0:  # Handle zero interest rate case
        return principal / term
    return (rate * principal) / (1 - (1 + rate) ** -term)

def calculate_principal_payment(mortgage_payment, mortgage_rate, outstanding_mortgage):
    """Calculate the principal portion of the mortgage payment."""
    interest_payment = outstanding_mortgage * mortgage_rate
    return mortgage_payment - interest_payment

# Bellman Equation Implementation
def optimize_allocation():
    optimal_allocations = np.zeros((YEARS,))
    
    for sim in range(SIMULATIONS):
        mortgage = START_MORTGAGE
        stocks = INITIAL_STOCKS
        equity = PROPERTY_VALUE - mortgage
        property_value = PROPERTY_VALUE
        
        earnings = INITIAL_EARNINGS
        spending = NON_MORTGAGE_SPENDING
        
        for year in range(YEARS):
            stock_return = stock_returns[sim, year]
            mortgage_rate = mortgage_rates[sim, year]

            mortgage_payment = calculate_pmt(mortgage_rate, MORTGAGE_TERM, mortgage)
            tax = tax_calculation(earnings)
            excess_earnings = calculate_excess_earnings(earnings, tax, spending, mortgage_payment)
            principal_payment = calculate_principal_payment(mortgage_payment, mortgage_rate, mortgage)

            # Try different allocations to find the optimal one
            best_allocation = 0
            best_utility = -np.inf
            for x in np.linspace(0, 1, 101):  # Test allocations from 0 to 1 in increments of 0.01
                stocks_next = stocks * (1 + stock_return) + x * excess_earnings
                mortgage_next = max(mortgage - principal_payment - (1 - x) * excess_earnings, 0)
                property_value_next = property_value * (1 + PROPERTY_GROWTH)

                equity_next = property_value_next - mortgage_next
                wealth_next = equity_next + stocks_next

                discounted_utility = utility(wealth_next, RISK_AVERSION) / ((1 + DISCOUNT_RATE) ** year)

                if discounted_utility > best_utility:
                    best_utility = discounted_utility
                    best_allocation = x

            # Update values
            stocks = stocks * (1 + stock_return) + best_allocation * excess_earnings
            mortgage = max(mortgage - principal_payment - (1 - best_allocation) * excess_earnings, 0)
            property_value = property_value * (1 + PROPERTY_GROWTH)
            
            optimal_allocations[year] += best_allocation
            
            # Update income and spending for next year
            earnings *= (1 + EARNINGS_GROWTH)
            spending *= (1 + SPENDING_GROWTH)

    # Average the optimal allocations across simulations
    optimal_allocations /= SIMULATIONS

    return optimal_allocations

# Run optimization
optimal_allocations = optimize_allocation()

# Plot the results
plt.figure(figsize=(10, 6))
plt.plot(range(1, YEARS + 1), optimal_allocations, marker='o')
plt.title("Optimal Stock Allocation Over Time")
plt.xlabel("Year")
plt.ylabel("Fraction of Excess Earnings Allocated to Stocks")
plt.grid()
plt.show()
```

## Answer by phdstudent (score 4)

https://quant.stackexchange.com/a/81518

It would be useful if you could show exactly the model you are trying to solve. In any case, solving a Bellman equation problem by monte carlo is a terrible idea. You would need a waaaay larger simulation.

You want to solve the belman equation first and then simulate the model. You need to approximate returns (look at Tauchen-Hussey approximation), and then solve the bellman by backward induction. Should be a trivial process I can sketch it but will need a better view of the problem you are trying to solve.