import numpy as np
def solve_hjb_trading(Q_total, T, dt, dq, eta, phi, alpha):
"""
Resolve numericamente a HJB para Execução Ótima / Ocultação de Ordem.
Parâmetros:
Q_total : Total de ações para liquidar
T : Tempo total de execução (ex: 1.0 = 1 dia)
dt : Passo de tempo
dq : Passo do inventário
eta : Parâmetro de impacto temporário (custo por velocidade)
phi : Penalidade de risco por inventário retido (urgência)
alpha : Penalidade por inventário não liquidado em T
"""
# Definição das grelhas temporais e de inventário
time_grid = np.arange(0, T + dt, dt)
inv_grid = np.arange(0, Q_total + dq, dq)
N_t = len(time_grid)
N_q = len(inv_grid)
# Matriz de Valor V[tempo, inventário]
V = np.zeros((N_t, N_q))
# Matriz de Política Ótima u[tempo, inventário] (velocidade de venda)
u_opt = np.zeros((N_t, N_q))
# Condição Terminal no tempo T (Penalidade se sobrar inventário)
V[-1, :] = -alpha * (inv_grid ** 2)
# Algoritmo de Indução Retrocedente (Backward Induction)
for t_idx in range(N_t - 2, -1, -1):
for q_idx in range(1, N_q):
q = inv_grid[q_idx]
# Aproximação de Diferença Finita para a derivada dV/dq (Backward Difference)
dv_dq = (V[t_idx + 1, q_idx] - V[t_idx + 1, q_idx - 1]) / dq
# Otimização HJB: u* = dv_dq / (2 * eta)
# Como estamos a vender, a velocidade reduz o inventário, forçamos u_star adequado
u_star = -dv_dq / (2 * eta)
# Limitar u_star para garantir estabilidade numérica e física
u_star = max(0, min(u_star, q / dt))
# Equação HJB Discretizada: V_t = V_next + dt * (maximizando custos)
# Custo de inventário retido (phi * q^2) reduz a utilidade
running_cost = - phi * (q ** 2) - eta * (u_star ** 2)
V[t_idx, q_idx] = V[t_idx + 1, q_idx] + dt * (running_cost - u_star * dv_dq)
u_opt[t_idx, q_idx] = u_star
return time_grid, inv_grid, u_opt, V
# --- Execução do Exemplo de Simulação ---
if __name__ == "__main__":
# Parâmetros de Mercado Fictícios
Q_init = 100000 # 100.000 ações
T_horizon = 1.0 # Janela de trading
dt_step = 0.01 # 100 passos de tempo
dq_step = 1000 # Resolução do inventário
eta_param = 0.05 # Custo do impacto de mercado imediato
phi_param = 0.1 # Aversão ao risco (evitar flutuações do mercado)
alpha_param = 2.0 # Penalidade severa para ordens não executadas no final
times, inventories, policy, value = solve_hjb_trading(
Q_init, T_horizon, dt_step, dq_step, eta_param, phi_param, alpha_param
)
print("📈 Tabela de Velocidade Ótima de Trading (Ações/dt) nos primeiros passos:")
print("Tempo \\ Inventário |", " | ".join([f"{int(q)} ac" for q in inventories[-4:]]))
print("-" * 65)
for t_i in range(5):
row = [f"{policy[t_i, q_i]:.2f}" for q_i in range(len(inventories)-4, len(inventories))]
print(f"t = {times[t_i]:.2f} | " + " | ".join(row))
