Here we provide another answer using Ito's calculus. It appears involved, but it also has its own interest.
Given the short rate dynamics
\begin{align*}
dr_t = \nu(r_t, t) dt + \rho(r_t, t) dW_t,
\end{align*}
we define the function
\begin{align*}
g(x, t, T) = -\ln E\left(e^{-\int_t^T r_s ds} \,\big|\, r_t = x\right).
\end{align*}
The forward rate $f(t, T)$ is then defined by
\begin{align*}
f(t, T) = \frac{\partial g}{\partial T}(r_t, t, T).
\end{align*}
Using Ito's lemma,
\begin{align*}
df(t, T) &= \frac{\partial^2 g}{\partial t \partial T} dt + \frac{\partial^2 g}{\partial r_t \partial T}dr_t + \frac{1}{2}\frac{\partial^3 g}{\partial^2 r_t \partial T}d\langle r, r\rangle_t\\
&=\left(\frac{\partial^2 g}{\partial t \partial T} + \frac{\partial^2 g}{\partial r_t \partial T} \nu(r_t, t) + \frac{1}{2}\rho(r_t, t)^2\frac{\partial^3 g}{\partial^2 r_t \partial T} \right)dt + \rho(r_t, t)\frac{\partial^2 g}{\partial r_t \partial T}dW_t\\
&= \sigma(t, T)\Sigma(t, T)dt + \sigma(t, T) dW_t,
\end{align*}
where
\begin{align*}
\sigma(t, T) &= \rho(r_t, t)\frac{\partial^2 g}{\partial r_t\partial T}(r_t, t, T), \\
\Sigma(t, T) &= \int_t^T\sigma(t, s) ds = \rho(r_t, t)\frac{\partial g}{\partial r_t}(r_t, t, T),\\
\sigma(t, T)\Sigma(t, T) &= \frac{\partial^2 g}{\partial t \partial T} + \frac{\partial^2 g}{\partial r_t \partial T} \nu(r_t, t) + \frac{1}{2}\rho(r_t, t)^2\frac{\partial^3 g}{\partial^2 r_t \partial T}.
\end{align*}
Then
\begin{align*}
\rho(r_t, t)^2\frac{\partial g}{\partial r_t}(r_t, t, T)\frac{\partial^2 g}{\partial r_t\partial T}(r_t, t, T) = \frac{\partial^2 g}{\partial t \partial T} + \frac{\partial^2 g}{\partial r_t \partial T} \nu(r_t, t) + \frac{1}{2}\rho(r_t, t)^2\frac{\partial^3 g}{\partial^2 r_t \partial T},
\end{align*}
that is,
\begin{align*}
\frac{1}{2}\rho(r_t, t)^2\frac{\partial }{\partial T}\left[\left(\frac{\partial g}{\partial r_t}\right)^2 \right] = \frac{\partial^2 g}{\partial t \partial T} + \frac{\partial^2 g}{\partial r_t \partial T} \nu(r_t, t) + \frac{1}{2}\rho(r_t, t)^2\frac{\partial^3 g}{\partial^2 r_t \partial T}.
\end{align*}
Note that
\begin{align*}
\int_t^T \frac{\partial }{\partial u}\left[\left(\frac{\partial g}{\partial r_t}(r_t, t, u)\right)^2 \right]du &= \left(\frac{\partial g}{\partial r_t}(r_t, t, T)\right)^2,\\
\int_t^T \frac{\partial^2 g}{\partial t \partial u}(r_t, t, u) du &= \lim_{s\rightarrow t+} \frac{\partial }{\partial t}\int_s^T \frac{\partial g}{\partial u}(r_t, t, u) du\\
&= \lim_{s\rightarrow t+} \frac{\partial }{\partial t} \big(g(r_t, t, T) -g(r_t,t, s)\big)\\
&=\frac{\partial g}{\partial t} +r_t.
\end{align*}
See Addendum below for more details. Then,
\begin{align*}
\frac{1}{2}\rho(r_t, t)^2\left(\frac{\partial g}{\partial r_t}\right)^2 = \frac{\partial g}{\partial t} +r_t + \frac{\partial g}{\partial r_t} \nu(r_t, t) + \frac{1}{2}\rho(r_t, t)^2\frac{\partial^2 g}{\partial^2 r_t}. \tag{1}
\end{align*}
Note that, we can also obtain Equation $(1)$ using the PDE for the bond price $ P$ (see PDE for Pricing Interest Rate Derivatives) and then make the substitution $P=e^{-g}$.
Moreover, note that
\begin{align*}
r_t &= f(t, t)\\
&=f(0, t) - \int_0^t \sigma(s, t)\Sigma(s, t)ds + \int_0^t \sigma(s, t) dW_s.
\end{align*}
In Ho-Lee model, $\rho(r_t, t) = \sigma$ and $\nu(r_t, t)=\theta_t$. Then for any $t>0$,
\begin{align*}
Var(r_t - f(t, t)) = \int_0^t E\left(\sigma(s, t)-\sigma \right)^2ds = 0
\end{align*}
That is,
\begin{align*}
\sigma(t, T) &= \sigma,\\
\frac{\partial^2 g}{\partial r_t\partial T}(r_t, t, T) &= 1, \\
\frac{\partial g}{\partial r_t}(r_t, t, T) &= T-t.
\end{align*}
From $(1)$,
\begin{align*}
\frac{1}{2} \sigma^2 (T-t)^2 = \frac{\partial g}{\partial t} +r_t + (T-t) \theta_t,
\end{align*}
Therefore,
\begin{align*}
\frac{\sigma^2}{6} (T-t)^3 = g(r_t, T, T) - g(r_t, t, T) +r_t(T-t) + \int_t^T (T-s) \theta_s ds.
\end{align*}
That is,
\begin{align*}
g(r_t, t, T) &= r_t(T-t) - \frac{\sigma^2}{6} (T-t)^3 + \int_t^T (T-s) \theta_s ds.
\end{align*}
Addendum
We note that
\begin{align*}
\lim_{s\rightarrow t+} \frac{\partial }{\partial t}g(r_t, t, s) &=\lim_{s\rightarrow t+}\lim_{\delta \rightarrow 0+}\frac{-\ln E\left(e^{-\int_{t+\delta}^s r_udu}\mid \mathcal{F}_{t+\delta} \right) + \ln E\left(e^{-\int_t^s r_udu}\mid \mathcal{F}_t \right)}{\delta}\\
&=\lim_{s\rightarrow t+}\lim_{\delta \rightarrow 0+}\frac{-\ln E\left(e^{\int_t^{t+\delta} r_udu}e^{-\int_t^s r_udu}\mid \mathcal{F}_{t+\delta} \right) + \ln E\left(e^{-\int_t^s r_udu}\mid \mathcal{F}_t \right)}{\delta}\\
&=\lim_{s\rightarrow t+}\lim_{\delta \rightarrow 0+}\frac{-\int_t^{t+\delta} r_udu -\ln E\left(e^{-\int_t^s r_udu}\mid \mathcal{F}_{t+\delta} \right) + \ln E\left(e^{-\int_t^s r_udu}\mid \mathcal{F}_t \right)}{\delta}\\
&=\lim_{\delta \rightarrow 0+}\lim_{s\rightarrow t+}\frac{-\int_t^{t+\delta} r_udu -\ln E\left(e^{-\int_t^s r_udu}\mid \mathcal{F}_{t+\delta} \right) + \ln E\left(e^{-\int_t^s r_udu}\mid \mathcal{F}_t \right)}{\delta}\\
&=-r_t.
\end{align*}