For the geometric Brownian motion $d X_t = \mu X_t dt + \sigma X_t d B_t$ being started at unity $X_0 = 1 $ the probability density function of the running maximum $M_T := sup(0<s<T: X_s)$ is given by the formula below:
\begin{eqnarray}
P_1 \left( M_T \in db \right)/d b &=&
b^{\frac{2 \mu }{\sigma ^2}-2} \left(\frac{\sqrt{\frac{2}{\pi }}
b^{\frac{1}{2}-\frac{\mu }{\sigma ^2}} \exp \left(-\frac{4
\log ^2(b)+T^2 \left(\sigma ^2-2 \mu \right)^2}{8 \sigma ^2
T}\right)}{\sigma \sqrt{T}}+\left(\frac{1}{2}-\frac{\mu
}{\sigma ^2}\right) \text{erfc}\left(\frac{\sqrt{\sigma ^2 T}
\left(\frac{\log (b)}{\sigma ^2 T}+\frac{\mu }{\sigma
^2}-\frac{1}{2}\right)}{\sqrt{2}}\right)\right)
\end{eqnarray}
where, following the notation from Ito&Mc Kean for example, $P_x \left(\omega\right) := P_x \left(\omega|X_0 = x\right)$.
Verification:
We start with the obvious identity $ P\left( M_T > b \right) = P \left( \tau_b < T \right ) $ where $\tau_b := inf(s<T: X_s=b)$ is the first hitting time of the barrier $b$.
The identity in question holds true as long as the process has no jumps.
Now take $x = 1$. We have :
\begin{eqnarray}
P_1 \left( M_T \in db \right)/d b &=&
-\frac{d}{d b} P_1 \left( \tau_b < T \right) \\
&=&
-\frac{d}{d b}
\int\limits_0^T
\frac{\left| \log \left(\frac{x}{b}\right)\right|
e^{-\frac{\left(\log \left(\frac{b}{x}\right)-\xi \left(\mu
-\frac{\sigma ^2}{2}\right)\right)^2}{2 \xi \sigma
^2}}}{\sqrt{2 \pi } \sqrt{\xi ^3} \sigma }
d\xi \\
&=&
-\frac{d}{d b}
\left[
\frac{1}{2} b^{\frac{\mu }{\sigma ^2}-\frac{1}{2}}
\left(
b^{-\left| \frac{\mu }{\sigma^2}-\frac{1}{2}\right| }
%
\text{erfc}\left(\frac{2 \log
(b)-T \left| \sigma ^2-2 \mu \right| }{2 \sqrt{2}
\sigma \sqrt{T}}\right)
%
+b^{\left| \frac{\mu }{\sigma
^2}-\frac{1}{2}\right| }
%
\text{erfc}\left(\frac{T
\left| \sigma ^2-2 \mu \right| +2 \log (b)}
{2\sqrt{2} \sigma \sqrt{T}}\right)
%
\right)
\right]
\end{eqnarray}
See Mathematica code snippet for numerical verification:
In[1124]:= Clear[T, s, b, a, \[Mu], \[Sigma], A, B, t];
(*This is the CDF of the first hitting time of barrier b when started \
from 1.*)
mcdf[b_] :=
1/2 b^(-(1/2) + \[Mu]/\[Sigma]^2) (
b^-Abs[-(1/
2) + \[Mu]/\[Sigma]^2] Erfc[(-T Abs[(-2 \[Mu] + \
\[Sigma]^2)] + 2 Log[b])/(2 Sqrt[2] \[Sigma]^1 Sqrt[T ])] +
b^+Abs[-(1/2) + \[Mu]/\[Sigma]^2] Erfc[(
T Abs[(-2 \[Mu] + \[Sigma]^2)] + 2 Log[b] )/(
2 Sqrt[2] \[Sigma]^1 Sqrt[T ])]);
{T, \[Mu], \[Sigma], b0} =
RandomReal[{1, 2}, 4, WorkingPrecision -> 50];
(*This is the CDF of the running supremum of a geometric Brownian \
motion with drift \[Mu] and variance (\[Sigma]^2).*)
b^(-2 + (2 [Mu])/[Sigma]^2) ( (1/2 - [Mu]/[Sigma]^2) Erfc[(
Sqrt[T [Sigma]^2] ( ([Mu]/[Sigma]^2 - 1/2) + Log[b]/(
T [Sigma]^2)))/ Sqrt[2] ] + (
b^(1/2 - [Mu]/[Sigma]^2)
E^(-((T^2 (-2 [Mu] + [Sigma]^2)^2 + 4 Log[b]^2)/(
8 T [Sigma]^2))) Sqrt[2/[Pi]])/(Sqrt[T] [Sigma])) /. b :> b0
-D[mcdf[b], b] /. b :> b0
Out[1127]= 0.259499378356891943946012789179083720641292239286
Out[1128]= 0.259499378356891943946012789179083720641292239286