I have the following problem: a standard Brownian motion Y(t) (starting at zero) is sampled at discrete times $t_j=jh$ (step size $h$). Given thresholds $d_j$ the process is said to die at time $t=t_j$ if $Y(t_j)<d_j$.
The thresholds $d_j$ are to be computed so that the probability that $Y(t)$ dies at time $t_j$ given that it is still alive at this time is a fixed probability $q$. The condition for the $d_j$ can then be stated as $$ P(Y(t_{j+1})<d_{j+1}\big|\inf\nolimits_{k\leq j}(Y(t_k)-d_k)\geq 0)=q $$ Using the Markov property this simplifies to $$ P(Y(t_{j+1})<d_{j+1}|Y(t_j)\geq d_j)=q $$ and suggests a recursive computation where we compute $d_{j+1}$ from $d_j$ by solving the following equation for $d=d_{j+1}$: $$ G(d,d_j)=q\quad\text{ where }\quad G(d,d_j)=P(Y(t_{j+1})<d\big|Y(t_j)\geq d_j)) $$ This function is increasing in the variable $d$ so the equation can be solved numerically by continued bisection once we have computed the function $G(d,d_j)$. To do this we use that the increment $Y(t_{j+1})-Y(t_j)$ is independent of $Y(t_j)$ and distributed as $N(0,h)$. Using also that $Y(t_j)\sim N(0,t_j)$ we can write
\begin{align*} G(d,d_j)&=P(Y(t_{j+1})-Y(t_j)<d-Y(t_j)\big|Y(t_j)\geq d_j)) \\&= \frac{1}{P(Y(t_j)\geq d_j)} \int_{d_j}^\infty P(Y(t_{j+1})-Y(t_j)<d-y\big|Y(t_j)=y)P_{Y(t_j)}(dy) \\&= \frac{1}{P(N(0,1)\geq d_j/\sqrt{t_j})} \int_{d_j}^\infty P(Y(t_{j+1})-y<d-y)P_{Y(t_j)}(dy) \\&= \frac{1}{1-F(d_j/\sqrt{t_j})} \int_{d_j}^\infty P\left(N(0,1)<\frac{d-y}{\sqrt h}\right)P_{Y(t_j)}(dy) \\&= \frac{1}{1-F(d_j/\sqrt{t_j})} \int_{d_j}^\infty F\left(\frac{d-y}{\sqrt h}\right)P_{Y(t_j)}(dy), \end{align*}
where $F$ is the standard normal cumulative distribution function. This integral is evaluated numerically. When the $d_j$ have the desired property then the probability that the process $Y(t)$ is alive at time $t=t_j$ is $(1-q)^j$.
However in subsequent simulation this is not the case. Instead $Y(t)$ dies more slowly than expected. Is there something wrong with the above reasoning?