Suppose $(V,g)$ is an $n$-dimensional oriented pseudo-inner product space. Let $\epsilon$ be the volume element. Then, the definition of the Hodge star is (up to sign, but this is the convention I’ll use and it gives the same answer as you write) defined such that for all $k$-forms $A,B$ on $V$,
\begin{align}
A\wedge \star B=\langle A,B\rangle_{g,k}\cdot\epsilon,
\end{align}
where $\langle\cdot,\cdot\rangle_{g,k}$ is the pseudo inner product induced on the space of $k$-forms $\bigwedge^k(V^*)$, and it is defined such that for all elementary wedges (i.e $a^1,\dots, a^k,b^1,\dots, b^k\in V^*$),
\begin{align}
\langle a^1\wedge\cdots\wedge a^k,b^1\wedge\cdots\wedge b^k\rangle_{g,k}&:=
\det\left(\langle a^i,b^j\rangle_{V^*}\right).
\end{align}
Ok so now we can get started with computations. Let’s say we have two arbitrary $2$-forms $A$ and $B$ on $V$. Fix a basis $\{v_1,\dots, v_n\}$ of $V$, and let $\{\phi^1,\dots, \phi^n\}$ be the dual basis for $V^*$. If we denote $g_{ij}=g(v_i,v_j)$, then the components of the inverse matrix are $g^{ij}=\langle\phi^i,\phi^j\rangle_{g,V^*}$ (here the inner product $\langle\cdot,\cdot\rangle_{g,V^*}$ is simply the one on $V^*$ induced by $g$).
Now, we can write
\begin{align}
A&=\frac{1}{2!}A_{ij}\,\phi^i\wedge\phi^j,\quad\text{and}\quad
B=\frac{1}{2!}B_{kl}\,\phi^k\wedge\phi^l.
\end{align}
So, using bilinearity of the inner product
\begin{align}
A\wedge \star B&=
\left\langle \frac{1}{2}A_{ij}\,\phi^i\wedge\phi^j,
\frac{1}{2}B_{kl}\,\phi^k\wedge\phi^l\right\rangle\cdot\epsilon\\
&=\frac{1}{4}A_{ij}B_{kl}\det
\begin{pmatrix}
\left\langle\phi^i,\phi^k\right\rangle_{g,V^*}&
\left\langle\phi^i,\phi^l\right\rangle_{g,V^*}\\
\left\langle\phi^j,\phi^k\right\rangle_{g,V^*}&
\left\langle\phi^j,\phi^l\right\rangle_{g,V^*}
\end{pmatrix}
\cdot\epsilon\\
&=\frac{1}{4}A_{ij}B_{kl}\det
\begin{pmatrix}
g^{ik}&g^{il}\\
g^{jk}&g^{jl}
\end{pmatrix}
\cdot\epsilon\\
&=\frac{1}{4}A_{ij}B_{kl}\cdot\left[g^{ik}g^{jl}-g^{il}g^{jk}\right]\cdot\epsilon\\
&=\left[\frac{1}{4}A_{ij}B^{ij}-\frac{1}{4}A_{ij}B^{ji}\right]\cdot\epsilon\\
&=\frac{1}{2}A_{ij}B^{ij}\cdot\epsilon,
\end{align}
where towards the end I did the usual index raising operations, and I used skew-symmetry of the components of $B$ in the last step.
This is almost the answer we seek. If we now apply the Hodge star to both sides of this equation, and keep in mind $\star\epsilon=(-1)^{\#}$, where $\#$ is the number of minus signs in the Sylvester form of the pseudo inner product $g$, then we get that
\begin{align}
\star(A\wedge \star B)&=\frac{(-1)^{\#}}{2}A_{ij}B^{ij}.
\end{align}
In the case of Lorentzian signature, we have $\#=1$ (assuming we take the one minus, mostly plus signature), so we get the $-\frac{1}{2}$ as expected. Now, you can specialize to $A=B$ if you want.
And you can easily generalize the argument to prove that for $k$-forms in general, $A\wedge\star B=\frac{1}{k!}A_IB^I\cdot\epsilon$, where the sum is over all $k$-tuples $I=(i_1,\dots, i_k)$ with each $i_{\alpha}\in\{1,\dots, n\}$; the $\frac{1}{k!}$ term accounts for the fact that we haven’t arranged for $I$ to be in increasing order.
