Dynamic Programming

from IPython.display import display, HTML
display(HTML("<style>.container { width:100% !important; }</style>"))

import numpy as np
import matplotlib.pyplot as plt
import math as m
from scipy import stats as st
from scipy import optimize
from scipy.optimize import minimize_scalar, bisect
from typing import NamedTuple
import time            # Imports system time module to time your script

plt.close('all')  # close all open figures

T = 9
beta = 0.96
kv = np.zeros(T+1,float)
cv = np.zeros(T+1,float)
uv = np.zeros(T+1,float)
kv[0] = 100  # k0
cv[0] = (1.0-beta)/(1.0-beta**(T+1)) * kv[0]  # c0
uv[0] = np.log(cv[0])

for i in range(1,T+1):
    #print "i=" + str(i)
    cv[i] = beta * cv[i-1]
    kv[i] = kv[i-1] - cv[i-1]

    # Period utility with discounting
    uv[i] = beta**i *np.log(cv[i])

np.sum(uv)  # total utility

print("cv = " + str(cv))
print("kv = " + str(kv))

fig, ax = plt.subplots(2,1)
plt.subplots_adjust(wspace=0.4, hspace=0.8)
#
ax[0].plot(cv, '-o')
ax[0].set_ylabel(r'$c^*_t$')
ax[0].set_xlabel('Period t')
ax[0].set_title('Optimal Consumption')
#
ax[1].plot(kv, '-o')
ax[1].set_ylabel(r'$k_t$')
ax[1].set_xlabel('Period t')
ax[1].set_title('Cake size')
#
plt.show()

def u(c, γ):

    return c**(1 - γ) / (1 - γ)
def v_star(k, β, γ):

    return (1 - β**(1 / γ))**(-γ) * u(k, γ)

β, γ = 0.95, 1.2
k_grid = np.linspace(0.1, 5, 100)

fig, ax = plt.subplots()

ax.plot(k_grid, v_star(k_grid, β, γ), label='value function')

ax.set_xlabel('$k$', fontsize=12)
ax.legend(fontsize=12)

plt.show()

import numpy as np
import matplotlib.pyplot as plt

# Example 1: T(x) = 0.5x
T = lambda x: 0.5 * x
x_vals = np.linspace(-10, 10, 200)

plt.figure(figsize=(6,4))
plt.plot(x_vals, T(x_vals), label="T(x) = 0.5x")
plt.plot(x_vals, x_vals, 'k--', label="45° line (fixed points)")
plt.title("Contraction Example 1: T(x) = 0.5x")
plt.xlabel("x")
plt.ylabel("T(x)")
plt.legend()
plt.grid(True)
plt.show()

# Example 2: T(x) = 0.5x + 2
T = lambda x: 0.5 * x + 2

# Iterate starting from any x0
x0 = 0
x_vals = [x0]
for _ in range(10):
    x_vals.append(T(x_vals[-1]))

print("Fixed point should be 4")
print("Iterates:", x_vals)

# VFI implementation for a simple consumption-savings problem with discrete income (runnable)
# This version has extra inline comments explaining shapes, broadcasting, and algorithmic steps.

import numpy as np
import matplotlib.pyplot as plt

# -------------------------
# Model parameters
# -------------------------
# Discount factor: how much agents value the future vs present
beta = 0.96

# Interest rate on assets (gross return = 1 + r)
r = 0.04

# CRRA risk aversion parameter (gamma = 1 is log utility; here we use 2 as a common example)
gamma = 2.0

# Asset grid: we discretize assets between a_min and a_max with Na points
a_min, a_max, Na = 0.0, 20.0, 200
a_grid = np.linspace(a_min, a_max, Na)  # shape (Na,)

# -------------------------
# Income process (discrete)
# -------------------------
# Two income states for simplicity: low and high
y_vals = np.array([0.5, 1.5])           # shape (Ny,) where Ny = 2

# Transition matrix P: P[iy, jy] = Prob(y' = jy | current y = iy)
# Rows should sum to 1. Here we use a simple persistence structure.
P = np.array([[0.9, 0.1],
              [0.1, 0.9]])              # shape (Ny, Ny)

# -------------------------
# Utility function
# -------------------------
def utility(c):
    """CRRA utility with safeguard to avoid taking powers/log of nonpositive numbers.
    This function accepts array input via numpy broadcasting and returns an array of same shape.
    """
    c_safe = np.maximum(c, 1e-12)  # prevent nonpositive consumption
    if gamma == 1.0:
        return np.log(c_safe)
    else:
        return (c_safe ** (1.0 - gamma)) / (1.0 - gamma)

# -------------------------
# Value function initialization
# -------------------------
# V[a_index, y_index] stores the value of being at assets a_grid[a_index] and income y_vals[y_index]
V = np.zeros((Na, len(y_vals)))  # initial guess: zeros

# Convergence settings
tol = 1e-6
max_iter = 2000

# -------------------------
# Value Function Iteration (VFI)
# -------------------------
# The main loop performs Bellman updates until V converges.
# At each iteration and for each current income state 'iy' we:
# 1. For each current asset a (index ia), consider all candidate next-period assets a' (indexed by ia_prime).
# 2. Compute feasible consumption c = (1+r)*a + y - a' for each pair (a, a').
#    This produces a matrix of shape (Na, Na) where rows correspond to current a and columns to a'.
# 3. Compute one-period utility U(c) for the entire (Na, Na) matrix.
# 4. Compute the expected continuation value for each candidate a' as:
#       E[V(a', y')] = sum_{y'} P[iy, y'] * V[a_prime_index, y']
#    This yields a vector of length Na (indexed by a').
# 5. Add the discounted continuation value to the one-period utility (broadcasted over rows).
# 6. Maximize over a' (columns) to update V_new[a, iy].
#
# Vectorization notes:
# - We build the consumption matrix c with broadcasting: a_grid[:, None] (shape (Na,1)) and a_grid[None, :] (shape (1,Na)).
# - The continuation value cont is computed as an expected value across next-period income states using matrix multiply.
# - cont has shape (Na,) and is added to every row of U (hence cont[None, :]).

for it in range(max_iter):
    V_new = np.empty_like(V)  # placeholder for updated value function

    # Loop over current income states (small dimension Ny = 2 here — looping is fine)
    for iy, y in enumerate(y_vals):
        # Build candidate consumption matrix:
        # For current assets along rows (shape (Na,1)) and a' along columns (shape (1,Na))
        # c[ia, ia_prime] = (1+r)*a_grid[ia] + y - a_grid[ia_prime]
        c = (1.0 + r) * a_grid[:, None] + y - a_grid[None, :]
        c = np.maximum(c, 1e-12)  # enforce nonnegativity for numerical stability

        # One-period utility for every (a, a') pair: shape (Na, Na)
        U = utility(c)

        # Expected continuation value for each candidate a' (vector of length Na):
        # P[iy, :] is the row of transition probabilities from current income state iy
        # V.T has shape (Ny, Na) after transpose, so (P[iy, :] @ V.T) yields a vector of length Na
        # Multiply by beta to discount future utility.
        cont = beta * (P[iy, :] @ V.T).T  # shape (Na,)

        # Right-hand side: one-period utility + discounted expected continuation value
        # We add cont[None, :] (shape (1, Na)) to U (shape (Na, Na)); broadcasting applies across rows.
        RHS = U + cont[None, :]

        # Maximize over columns (a' choices) for each row (current a)
        # axis=1 corresponds to maximizing across candidate a' for each current a
        V_new[:, iy] = np.max(RHS, axis=1)

    # Compute sup-norm difference for convergence check
    diff = np.max(np.abs(V_new - V))
    V[:] = V_new  # update value function

    if diff < tol:
        print(f'Converged after {it+1} iterations, diff={diff:.2e}')
        break
else:
    print('VFI did not converge within max_iter')

# -------------------------
# Recovering the policy function
# -------------------------
# After convergence, recover the argmax (policy) that attains the maximum for each (a, y).
# policy_idx[a_index, y_index] = index of a' chosen (an integer index into a_grid)
policy_idx = np.empty_like(V, dtype=int)

for iy, y in enumerate(y_vals):
    # Recompute the same RHS components to obtain argmax (this repeats some computation;
    # in production code you might store argmax during the VFI loop to avoid recomputing)
    c = (1.0 + r) * a_grid[:, None] + y - a_grid[None, :]
    c = np.maximum(c, 1e-12)
    U = utility(c)
    cont = beta * (P[iy, :] @ V.T).T
    RHS = U + cont[None, :]
    policy_idx[:, iy] = np.argmax(RHS, axis=1)   # index of optimal a' for each current a

# Convert index to actual asset choices a' (policy in levels)
policy_a_prime = a_grid[policy_idx]  # shape (Na, Ny)

# -------------------------
# Plots: value functions and policy functions
# -------------------------
fig, axes = plt.subplots(1, 2, figsize=(10, 4))

# Value functions for the two income states
axes[0].plot(a_grid, V[:, 0], label=f'y={y_vals[0]}')
axes[0].plot(a_grid, V[:, 1], label=f'y={y_vals[1]}')
axes[0].set_title('Value function V(a,y)')
axes[0].set_xlabel('Assets a')
axes[0].set_ylabel('Value V')
axes[0].legend()

# Policy: a' choice as function of current assets for each income state
axes[1].plot(a_grid, policy_a_prime[:, 0], label=f'y={y_vals[0]}')
axes[1].plot(a_grid, policy_a_prime[:, 1], label=f'y={y_vals[1]}')
axes[1].set_title("Policy: a' (next assets)")
axes[1].set_xlabel('Assets a')
axes[1].set_ylabel("Chosen a'")
axes[1].legend()

plt.tight_layout()
plt.show()

# EGM vs VFI: implementations and timing comparison
import numpy as np
import time
import matplotlib.pyplot as plt
from scipy import interpolate

# -------------------------
# Model primitives (deterministic income)
# -------------------------
beta = 0.96       # discount factor
r = 0.04          # interest rate
gamma = 2.0       # CRRA risk aversion (gamma=1 is log utility)
y = 1.0           # deterministic income each period
a_min = 0.0       # borrowing constraint (natural boundary)
# We'll vary a_max and Na (grid size) below for timing comparison.

# Utility and marginal utility functions
def u(c):
    c_safe = np.maximum(c, 1e-12)
    if gamma == 1.0:
        return np.log(c_safe)
    else:
        return (c_safe ** (1.0 - gamma)) / (1.0 - gamma)

def mu(c):
    # marginal utility u'(c) = c^{-gamma} for CRRA
    c_safe = np.maximum(c, 1e-12)
    return c_safe ** (-gamma)

def mu_inv(mu_val):
    # inverse marginal utility: given marginal utility, return c
    # For CRRA: mu = c^{-gamma} => c = mu^{-1/gamma}
    return mu_val ** (-1.0 / gamma)

# -------------------------
# VFI implementation for deterministic model (function)
# -------------------------
def solve_vfi_det(a_grid, beta=beta, r=r, y=y, tol=1e-6, max_iter=2000):
    """Solve deterministic consumption-savings via naive VFI on given asset grid.
    Returns: consumption policy c_policy (shape Na,), value function V (Na,), and timing info.
    """
    Na = len(a_grid)
    V = np.zeros(Na)                # initial guess for value function
    c_policy = np.zeros(Na)         # to store policy in levels
    # Precompute resource (cash-on-hand) for state a and candidate a'
    # We'll vectorize across a (rows) and a' (cols)
    A = a_grid[:, None]             # shape (Na,1)
    Aprime = a_grid[None, :]        # shape (1,Na)

    start = time.perf_counter()
    for it in range(max_iter):
        V_new = np.empty_like(V)
        # compute RHS for all (a, a') pairs
        cons = (1.0 + r) * A + y - Aprime   # shape (Na, Na)
        cons = np.maximum(cons, 1e-12)
        U = u(cons)                         # shape (Na, Na)
        # continuation value is V evaluated at a' (vector length Na), broadcast across rows
        cont = beta * V[None, :]            # shape (1, Na)
        RHS = U + cont                      # shape (Na, Na)
        V_new = np.max(RHS, axis=1)        # max over columns (a') for each row (current a)
        # Check convergence
        diff = np.max(np.abs(V_new - V))
        V[:] = V_new
        if diff < tol:
            break
    end = time.perf_counter()
    # Recover policy (argmax)
    cons = (1.0 + r) * A + y - Aprime
    cons = np.maximum(cons, 1e-12)
    U = u(cons)
    cont = beta * V[None, :]
    RHS = U + cont
    idx = np.argmax(RHS, axis=1)
    c_policy = cons[np.arange(Na), idx]
    return c_policy, V, end - start, it+1

# -------------------------
# EGM implementation for deterministic model (function)
# -------------------------
def solve_egm_det(a_grid, beta=beta, r=r, y=y, tol=1e-6, max_iter=2000):
    """Solve deterministic consumption-savings via EGM with safe interpolation."""
    Na = len(a_grid)
    # initial guess: consume half of cash-on-hand
    c_policy = 0.5 * ((1.0 + r) * a_grid + y)
    c_policy = np.maximum(c_policy, 1e-12)

    start = time.perf_counter()
    for it in range(max_iter):
        # 1) next-period consumption on a' grid
        c_next = c_policy.copy()
        # 2) next-period marginal utility
        muc_next = mu(c_next)
        # 3) invert Euler: u'(c) = beta(1+r)u'(c')
        factor = beta * (1.0 + r) * muc_next
        # guard against numerical under/overflow
        factor = np.maximum(factor, 1e-300)
        c_current_at_a_prime = mu_inv(factor)

        # 4) implied current asset that leads to (c_current, a')
        a_implied = (c_current_at_a_prime + a_grid - y) / (1.0 + r)

        # 5) borrowing constraint handling: if a_implied < a_min, the correct solution is a=a_min
        mask_below = a_implied < a_min
        if np.any(mask_below):
            a_implied[mask_below] = a_min
            # at the constraint the optimal a' is a_min; consumption there is c = (1+r)a_min + y - a'
            # since we are sitting at a=a_min, set a' = a_min for the constraint point
            # i.e., c = y + r*a_min - a_min = y + (r-1)*a_min; with a_min=0 this is y
            c_current_at_a_prime[mask_below] = (1.0 + r) * a_min + y - a_min

        # 6) sort by a_implied and collapse duplicate x's (take max c for each unique a)
        sort_idx = np.argsort(a_implied)
        x = a_implied[sort_idx]
        y_c = c_current_at_a_prime[sort_idx]

        # collapse duplicates robustly
        unique_x = []
        unique_c = []
        i = 0
        n = x.size
        while i < n:
            j = i + 1
            while j < n and x[j] == x[i]:
                j += 1
            unique_x.append(x[i])
            # at duplicated a (notably a_min), keep the largest feasible consumption
            unique_c.append(np.max(y_c[i:j]))
            i = j

        x = np.asarray(unique_x)
        y_c = np.asarray(unique_c)

        # ensure we have the constraint anchor explicitly (use correct c at a_min)
        if x[0] > a_min:
            x = np.insert(x, 0, a_min)
            y_c = np.insert(y_c, 0, (1.0 + r) * a_min + y - a_min)

        # if any tiny non-monotonicity remains, enforce strictly increasing x by jittering eps
        # (rarely needed, but keeps interp1d happy)
        eps = 0.0
        for k in range(1, x.size):
            if x[k] <= x[k-1]:
                eps = max(eps, 1e-12 * (1.0 + abs(x[k-1])))
                x[k] = x[k-1] + eps

        # 7) interpolate back to the original grid
        interp = interpolate.interp1d(x, y_c, kind='linear', bounds_error=False, fill_value='extrapolate')
        c_new = interp(a_grid)
        c_new = np.maximum(c_new, 1e-12)

        # 8) convergence
        diff = np.max(np.abs(c_new - c_policy))
        c_policy[:] = c_new
        if diff < tol:
            break

    end = time.perf_counter()
    return c_policy, end - start, it + 1

# -------------------------
# Timing comparison across grid sizes
# -------------------------
grid_sizes = [100, 200, 400, 800]  # adjust depending on runtime tolerance in this environment
vfi_times = []
egm_times = []
vfi_iters = []
egm_iters = []

for Na in grid_sizes:
    a_max = 50.0
    a_grid = np.linspace(a_min, a_max, Na)

    # Time VFI (deterministic)
    c_vfi, V_vfi, t_vfi, it_vfi = solve_vfi_det(a_grid)
    vfi_times.append(t_vfi)
    vfi_iters.append(it_vfi)

    # Time EGM
    c_egm, t_egm, it_egm = solve_egm_det(a_grid)
    egm_times.append(t_egm)
    egm_iters.append(it_egm)

    print(f"Na={Na}: VFI time={t_vfi:.4f}s (iters={it_vfi}), EGM time={t_egm:.4f}s (iters={it_egm})")

# -------------------------
# Plot timing results
# -------------------------
plt.figure(figsize=(6,4))
plt.plot(grid_sizes, vfi_times, marker='o', label='VFI (deterministic)')
plt.plot(grid_sizes, egm_times, marker='o', label='EGM')
plt.xlabel('Number of asset grid points (Na)')
plt.ylabel('Compute time (seconds)')
plt.title('Runtime comparison: VFI vs EGM (deterministic model)')
plt.legend()
plt.grid(True)
plt.show()

# -------------------------
# Show policy functions for the finest grid
# -------------------------
Na = grid_sizes[-1]
a_grid = np.linspace(a_min, 50.0, Na)
c_vfi, V_vfi, t_vfi, it_vfi = solve_vfi_det(a_grid)
c_egm, t_egm, it_egm = solve_egm_det(a_grid)

plt.figure(figsize=(8,4))
plt.plot(a_grid, c_vfi, label='VFI consumption policy')
plt.plot(a_grid, c_egm, label='EGM consumption policy', linestyle='--')
plt.xlabel('Assets a')
plt.ylabel('Consumption c(a)')
plt.title(f'Policy functions (Na={Na}) — VFI vs EGM')
plt.legend()
plt.grid(True)
plt.show()

from IPython.display import HTML, display

display(HTML("""
<style>
/* Expand notebook width */
.container { 
    width: 100% !important;
}

/* Disable output scrolling */
.output_scroll {
    height: auto !important;
    max-height: none !important;
    overflow-y: visible !important;
}

/* Also expand output area width */
.output_wrapper, .output {
    overflow-x: visible !important;
}
</style>
"""))

# EGM for stochastic income (discrete Markov) + demo run
import numpy as np
import time
import matplotlib.pyplot as plt
from scipy import interpolate

# Model primitives (shared)
beta = 0.96
r = 0.04
gamma = 2.0
a_min = 0.0

def mu(c):
    c_safe = np.maximum(c, 1e-12)
    return c_safe ** (-gamma)

def mu_inv(mu_val):
    return mu_val ** (-1.0 / gamma)

# -------------------------
# Stochastic EGM solver
# -------------------------
def solve_egm_stochastic(a_grid, y_vals, P, beta=beta, r=r, tol=1e-6, max_iter=2000):
    """Endogenous Grid Method for a model with discrete Markov income.
    Inputs:
      - a_grid: 1D array of asset grid points (Na,)
      - y_vals: 1D array of income states (Ny,)
      - P: transition matrix shape (Ny, Ny) with rows summing to 1
    Returns:
      - c_policy: consumption policy array shape (Na, Ny) mapping (a, y) -> c
      - elapsed_time: seconds
      - n_iter: iterations taken
    Notes:
      - We iterate on the consumption policy c(a, y).
      - For each current income state i and for each a' in the grid we compute expected
        marginal utility using P[i, :], invert Euler, get implied a, then interpolate.
    """
    Na = len(a_grid)
    Ny = len(y_vals)

    # Initial guess: consume a fraction of cash-on-hand for each income state
    # cash on hand = (1+r)*a + y, choose to consume half of it initially
    coh = (1.0 + r) * a_grid[:, None] + y_vals[None, :]   # shape (Na, Ny)
    c_policy = 0.5 * coh                                   # initial c(a,y), shape (Na, Ny)
    c_policy = np.maximum(c_policy, 1e-12)

    start = time.perf_counter()
    for it in range(max_iter):
        c_new = np.empty_like(c_policy)
        # For each current income state y_i compute implied a and c from a' grid
        for iy, y in enumerate(y_vals):
            # For each candidate a' (index j), we need c_next at (a'=a_grid[j], every y')
            # Extract c_next(a', y') as a vector of length Ny for each a' -- that's c_policy[j, :].
            # We'll loop over a' index vectorized with array ops.
            # c_next_matrix has shape (Na, Ny) where row j corresponds to a'=a_grid[j] across y'
            c_next_matrix = c_policy.copy()  # shape (Na, Ny): rows are a' indices

            # Compute marginal utility at next period for each (a', y') -> shape (Na, Ny)
            muc_next = mu(c_next_matrix)     # shape (Na, Ny)

            # Compute expected marginal utility conditional on current income state iy:
            # For each a' (row), take inner product with P[iy, :] over columns (y' dimension)
            # This yields a vector of length Na: muc_exp[j] = sum_{y'} P[iy, y'] * muc_next[j, y']
            muc_exp = muc_next @ P[iy, :].T   # shape (Na,)   (matrix-vector multiplication)

            # Invert Euler condition for each a' to get current consumption associated with that a'
            factor = beta * (1.0 + r) * muc_exp
            c_current_at_aprime = mu_inv(factor)   # shape (Na,)

            # Get implied current asset for each a' using budget identity:
            # a = (c + a' - y) / (1+r)
            a_implied = (c_current_at_aprime + a_grid - y) / (1.0 + r)  # shape (Na,)

            # Handle borrowing constraint: where a_implied < a_min, set a_implied = a_min and recompute c
            mask = a_implied < a_min
            if np.any(mask):
                a_implied[mask] = a_min
                # recompute c for these a' from budget: c = (1+r)*a_min + y - a'
                c_current_at_aprime[mask] = (1.0 + r) * a_min + y - a_grid[mask]

            # Now we have pairs (a_implied, c_current_at_aprime) of length Na.
            # Sort by a_implied (necessary for monotone interpolation)
            sort_idx = np.argsort(a_implied)
            a_sorted = a_implied[sort_idx]
            c_sorted = c_current_at_aprime[sort_idx]

            # Remove duplicate a's to avoid zero denominators in slope calculation
            a_unique, unique_idx = np.unique(a_sorted, return_index=True)
            c_unique = c_sorted[unique_idx]

            # Interpolate onto original a_grid to get c(a_grid, y_i)
            interp = interpolate.interp1d(a_unique, c_unique, kind='linear',
                                          bounds_error=False, fill_value='extrapolate')
            c_on_grid = interp(a_grid)
            c_on_grid = np.maximum(c_on_grid, 1e-12)
            c_new[:, iy] = c_on_grid

        # Convergence check across all (a,y)
        diff = np.max(np.abs(c_new - c_policy))
        c_policy[:] = c_new
        if diff < tol:
            break

    elapsed = time.perf_counter() - start
    return c_policy, elapsed, it+1

# -------------------------
# Demo run: small Markov income process
# -------------------------
y_vals = np.array([0.5, 1.5])   # low and high income
P = np.array([[0.9, 0.1],
              [0.1, 0.9]])      # transition matrix (rows sum to 1)

Na = 400
a_max = 50.0
a_grid = np.linspace(a_min, a_max, Na)

c_stoch_egm, t_stoch_egm, it_stoch_egm = solve_egm_stochastic(a_grid, y_vals, P)
print(f'Stochastic EGM converged in {it_stoch_egm} iters, time {t_stoch_egm:.4f}s')

# Plot policy functions for each income state
plt.figure(figsize=(8,4))
plt.plot(a_grid, c_stoch_egm[:,0], label=f'EGM c(a, y_low={y_vals[0]})')
plt.plot(a_grid, c_stoch_egm[:,1], label=f'EGM c(a, y_high={y_vals[1]})', linestyle='--')
plt.xlabel('Assets a')
plt.ylabel('Consumption c(a,y)')
plt.title('Stochastic EGM: Consumption policy by income state')
plt.legend()
plt.grid(True)
plt.show()

## 8. (Optional) Minimal Python illustration: value iteration for a small MDP

# simple value iteration for a small finite MDP
import numpy as np

# example: 2 states, 2 actions
P = {
  0: {0: [(1.0,0)], 1: [(1.0,1)]},   # action 0 keeps state 0, action 1 moves to state 1
  1: {0: [(1.0,1)], 1: [(1.0,1)]}    # absorbing state 1
}
R = {
  (0,0): 1.0, (0,1): 2.0,
  (1,0): 0.0, (1,1): 0.0
}
beta = 0.9
V = np.zeros(2)

for it in range(100):
    V_new = V.copy()
    for s in [0,1]:
        vals = []
        for a in [0,1]:
            expV = sum(p * V[s2] for (p,s2) in P[s][a])
            vals.append(R[(s,a)] + beta*expV)
        V_new[s] = max(vals)
    if np.max(np.abs(V_new - V)) < 1e-8:
        break
    V = V_new
print("Value function:", V)