Task: self_exciting_pricing_dp
Category: Operations research
Difficulty: 2
Tags: dynamic programming,finite difference,layer-parallel,pricing
Source: Yuan, Du \
Description: Hu 2026, "Dynamic Pricing Under Self-Exciting Arrival Processes," Operations Research: https://doi.org/10.1287/opre.2024.1172 & Finite-difference dynamic-pricing DP under self-exciting point-process demand.

Prompt template (YAML):
task_description: |
  The problem is to compute the optimal dynamic pricing policy for a product under a self-exciting demand arrival process (Hawkes process). 
  The problem is formulated as a continuous-time, continuous-state dynamic program (HJB equation). 
  We discretize the state space (excitement level x) and time to solve the DP backward in time. 
  At each time step, for each state, the algorithm must find the optimal purchase probability that maximizes the sum of immediate profit and the expected future profit from the jump in excitement level.

interface: |
  void solution_compute(
      int N_x, int N_T, int N_lambda,
      float dx, float dt, float alpha, float beta, float phi,
      float c, float a, float gamma, float k,
      float lambda_min, float lambda_max,
      float* V_out
  );

input_size_notes: |
  - N_x: Number of states (up to 100,000)
  - N_T: Number of time steps (up to 1,000)
  - N_lambda: Number of lambda values to search (up to 200)

output_contract: |
  - `V_out` is a pre-allocated host array of size `N_x`.
  - The function `solution_compute` must compute the final value function `V_prev` at t=0 and store it in `V_out`.

algorithm_background: |
  1. Initialize the value function array `V_prev` of size N_x with zeros.
  2. Loop backward in time for N_T steps:
  3. For each state index i in [0, N_x - 1] (in parallel):
  4.   Calculate x = i * dx.
  5.   Calculate the spatial derivative: V_deriv = (i == 0) ? 0.0 : (V_prev[i] - V_prev[i-1]) / dx.
  6.   Calculate the jump index: idx_alpha = min(i + round(alpha / dx), N_x - 1).
  7.   Calculate the jump value: V_jump = V_prev[idx_alpha] - V_prev[i].
  8.   Initialize max_val = -infinity (e.g., -1e30f).
  9.   For each lambda in the discretized set Lambda (N_lambda values linearly spaced from lambda_min to lambda_max):
  10.    Calculate val = r(lambda) + lambda * V_jump, where r(lambda) = lambda * (log(a / lambda) - c).
  11.    If val > max_val, update max_val = val.
  12.  Update the value function: V_new[i] = V_prev[i] + dt * (-x * beta * V_deriv + (k * x^gamma + phi) * max_val).
  13. Swap V_prev and V_new pointers.
  14. After the time loop, copy the final V_prev array to V_out.

hints_l2: |
  - You can precompute the lambda values and r(lambda) at the start of `solution_compute` and store them in device memory.
  - The outer loop over time steps is sequential, but the inner loop over states `i` can be parallelized across GPU threads.
  - Use double buffering for `V_prev` and `V_new` on the device to avoid allocations inside the time loop.

hints_l1: |
  - Allocate device memory for `V_prev` and `V_new` at the start of `solution_compute`.
  - In `solution_compute`, launch a kernel for each time step `t`.
  - In the kernel, each thread computes the update for one state `i`.
  - Since `V_prev` is read by multiple threads, consider using `__ldg()` or shared memory if appropriate, though L1 cache usually handles this well.
  - The lambda search loop is small (up to 200 iterations) and can be done sequentially by each thread.


CPU reference implementation (C):
#include <stdlib.h>
#include <stdio.h>
#include <math.h>

#ifdef __cplusplus
extern "C" {
#endif

void solution_compute(
    int N_x, int N_T, int N_lambda,
    float dx, float dt, float alpha, float beta, float phi,
    float c, float a, float gamma, float k,
    float lambda_min, float lambda_max,
    float* V_out
) {
    float* V_prev = (float*)calloc(N_x, sizeof(float));
    float* V_new = (float*)calloc(N_x, sizeof(float));
    
    float* lambdas = (float*)malloc(N_lambda * sizeof(float));
    float* r_lambdas = (float*)malloc(N_lambda * sizeof(float));
    
    float lambda_step = (lambda_max - lambda_min) / (N_lambda - 1);
    for (int j = 0; j < N_lambda; j++) {
        lambdas[j] = lambda_min + j * lambda_step;
        r_lambdas[j] = lambdas[j] * (logf(a / lambdas[j]) - c);
    }
    
    int alpha_idx_shift = (int)roundf(alpha / dx);
    
    for (int t = 0; t < N_T; t++) {
        for (int i = 0; i < N_x; i++) {
            float x = i * dx;
            float V_deriv = (i == 0) ? 0.0f : (V_prev[i] - V_prev[i-1]) / dx;
            
            int idx_alpha = i + alpha_idx_shift;
            if (idx_alpha > N_x - 1) idx_alpha = N_x - 1;
            
            float V_jump = V_prev[idx_alpha] - V_prev[i];
            
            float max_val = -1e30f;
            for (int j = 0; j < N_lambda; j++) {
                float val = r_lambdas[j] + lambdas[j] * V_jump;
                if (val > max_val) {
                    max_val = val;
                }
            }
            
            float h_x = k * powf(x, gamma);
            
            V_new[i] = V_prev[i] + dt * (-x * beta * V_deriv + (h_x + phi) * max_val);
        }
        
        // Swap
        float* temp = V_prev;
        V_prev = V_new;
        V_new = temp;
    }
    
    for (int i = 0; i < N_x; i++) {
        V_out[i] = V_prev[i];
    }
    
    free(V_prev);
    free(V_new);
    free(lambdas);
    free(r_lambdas);
}

void solution_free(void) {
    // No persistent state to free
}

#ifdef __cplusplus
}
#endif