I have proved in what follows that $\ast\ast\alpha=(-1)^{kn+k}\alpha$ for every $k$-form $\alpha$ on $\mathbb{R}^n$ and I would appreciate if someone would check my proof and/or point out how to improve it, thanks.
Proof. Let $\alpha$ be a $k$-form on $\mathbb{R}^n$, $\alpha=\sum\limits_{I}f_I dx_I$; we recall that the Hodge star operator is such that $\ast\alpha=\sum\limits_{I}f_I\varepsilon_{I}dx_{I^c}$, where $I^c$ denotes the complementary increasing multi-index, which consists of all the numbers between $1$ and $n$ that do not occur in $I$ and $\varepsilon_I=\begin{cases}1 & \text{if } dx_{I}dx_{I^c}=dx_1dx_2\cdots dx_n,\\ -1 & \text{if }dx_{I}dx_{I^c}=-dx_1dx_2\cdots dx_n \end{cases}$.
There are two cases two consider: $dx_{I}dx_{I^c}=dx_1dx_2\cdots dx_n$ and $dx_{I}dx_{I^c}=-dx_1dx_2\cdots dx_n.$
Suppose that $dx_{I}dx_{I^c}=dx_1dx_2\cdots dx_n$: then $\ast\alpha=\sum\limits_{I}f_I dx_{I^c}$ and $\ast\ast\alpha=\sum\limits_{I}f_I \varepsilon_{I^c}dx_{I}$.
Now, \begin{align}dx_{I^c}dx_{I}&\overset{(1)}{=}(-1)^{(n-k)k}dx_{I}dx_{I^c}\overset{(2)}{=}(-1)^{kn-k^2}dx_1dx_2\cdots dx_n\\ &\overset{(3)}{=}(-1)^{kn+k}dx_1dx_2\cdots dx_n\end{align} so $\varepsilon_{I^c}=(-1)^{kn+k}$ thus $**\alpha=\sum\limits_{I}f_I(-1)^{kn+k}dx_I=(-1)^{kn+k}\sum\limits_{I}f_Idx_I=(-1)^{kn+k}\alpha$, as desired.
Suppose instead that $dx_{I}dx_{I^c}=-dx_1dx_2\cdots dx_n$: then $\ast\alpha=\sum\limits_{I}f_I (-1) dx_{I^c},\ \ast\ast\alpha=\sum\limits_{I}f_I (-1)\varepsilon_{I^c}dx_{I}$ and since \begin{align}dx_{I^c}dx_{I}&\overset{(1)}{=}(-1)^{(n-k)k}dx_{I}dx_{I^c}\overset{(4)}{=}(-1)^{kn-k^2}(-1)dx_{1}dx_{2}\cdots dx_{n}\overset{(3)}{=}(-1)^{kn+k}(-1)dx_{1}dx_{2}\cdots dx_n\\&=(-1)^{kn+k+1}dx_1dx_2\cdots dx_{n}\end{align} it follows that $\varepsilon_{I^c}=(-1)^{kn+k+1}$ thus
\begin{align}\ast\ast\alpha&=\sum\limits_{I}f_I (-1)\varepsilon_{I^c}dx_{I}=\sum\limits_{I}f_I(-1)(-1)^{kn+k+1}dx_I=\sum\limits_{I}f_I(-1)^2(-1)^{kn+k}dx_I\\ &=(-1)^{kn+k}\sum\limits_{I}f_Idx_I=(-1)^{kn+k}\alpha,\end{align} as desired. $\square$
(1) we are shifting each of the $n-k\ \ dx_j$s in $dx_{I^c}$ through the $k\ \ dx_i$s in $dx_{I}$;
(2) multiplying out the exponent and from our hypothesis that $dx_{I}dx_{I^c}=dx_1dx_2\cdots dx_n$;
(3) $(-1)^{kn-k^2}=(-1)^{kn}(-1)^{-k^2}=(-1)^{kn}\frac{1}{(-1)^{k^2}}$ and since $k^2$ is even/odd if and only if $k$ is even/odd we get $\frac{1}{(-1)^{k^2}}=\frac{1}{(-1)^k}$ and being $\frac{1}{(-1)^k}=(-1)^k$ for all $k\in\mathbb{N}$ this equals $(-1)^{kn+k}$;
(4) multiplying out the exponent and from our hypothesis that $dx_{I}dx_{I^c}=-dx_1dx_2\cdots dx_n$;
Seems like it's on the right track. But one way to streamline the presentation, and avoid the case work, is to write the sign as $\epsilon_{I,I^c}=\text{sgn}(I,I^c)$, the sign of the permutation $(1,\dots, n)\mapsto (i_1,\dots, i_k,i_1',\dots, i_{n-k}')$, where $I=(i_1,\dots, i_k)$ and $I^c=(i_1',\dots, i_{n-k}')$ in strict increasing order. Then, we have \begin{align} \star\star(dx_I)=\star(\epsilon_{I,I^c}dx_{I^c})=\epsilon_{I,I^c}\star(dx_{I^c})&= \epsilon_{I,I^c}\epsilon_{I^c,(I^c)^c}\,dx_{(I^c)^c}\\ &=\epsilon_{I,I^c}\epsilon_{I^c,I}\,dx_I\\ &=(-1)^{k(n-k)}[\epsilon_{I,I^c}]^2\,dx_I\\ &=(-1)^{k(n-k)}\,dx_I. \end{align} So really I'm not doing anything different; the factor $(-1)^{k(n-k)}$ comes from making that many swaps to change $\epsilon_{I^c,I}$ into $\epsilon_{I,I^c}$, but this just avoids any case distinction. Then, by linearity it follows for all $k$-forms.