The answers to these questions are known, but, perhaps, not well-known. The typical approach is to first divide only by the local diffeomorphisms $\phi:\mathbb{R}^n\to\mathbb{R}^n$ that fix the origin and for which $\phi'(0):\mathbb{R}^n\to\mathbb{R}^n$ is the identity. This quotient is essentially sectioned by geodesic normal coordinates $(x^i)$, in which a given metric $g=g_{ij}(x)\,\mathrm{d}x^i\mathrm{d}x^j$ satisfies
$$
g_{ij}(x)\,x^j = g_{ij}(0)\,x^j.
$$
Such a $g$ can be expanded in Taylor series
$$
g_{ij}(x) = g_{ij}(0) + g^1_{ij}(x) + \cdots + g^k_{ij}(x) + \cdots
$$
where $g^k_{ij}=g^k_{ji}$ is homogeneous in $(x^i)$ of degree $k$. For $k>0$, the dimension of the vector space of such $\bigl(g^k_{ij}(x)\bigr)$ that satisfy the relations $g^k_{ij}(x)\,x^j=0$ is easily computed to be
$$
D_{k}(n) = {n\choose2}\,{{n+k-1}\choose{k}}-n\,{{n+k-1}\choose{k+1}}=\frac{n\,(k{-}1)}2\,{{n+k-1}\choose{k+1}},
$$
while $D_0(n) = {{n+1}\choose2}$. The group $\mathrm{GL}(n,\mathbb{R})$, as linear changes of the coordinates $(x^i)$, acts on the sum of these vector spaces. Since $g_{ij}(0)$ is positive definite, the action is transitive on the $0$-degree piece, with stabilizer $\mathrm{O}(n)\subset\mathrm{GL}(n,\mathbb{R})$, so this piece can be discarded and one can assume that $g_{ij}(0) = \delta_{ij}$.
Thus, one can regard the quotient space of diffeomorphism-equivalence classes of $k$-jets for $k\ge 2$ as the quotient of the sum of the $j$-th homogeneous pieces for $2\le j\le k$ divided by a natural action of $\mathrm{O}(n)$. For $n>2$, the generic element has a finite stabilizer, so the dimension of the named quotient space is
$$
\dim (S^{(k)}/\mathrm{Diff}_{k+1}) = D_2(n) + \cdots + D_k(n) - \dim\mathrm{O}(n),
$$
for $k\ge 2$.
In particular, the Poincaré series of the graded vector space of the positive degree terms, i.e., $D_2(n)\,t^2 +D_3(n)\,t^3 + \cdots$ is a rational function $p_n(t)/(1{-}t)^n$,
where $p_n(t)$ is a polynomial of degree $n$ with leading term $(-1)^n{n\choose2}\,t^n$, corresponding to the fact that, in Cartan's sense, Riemannian metrics in dimension $n$ depend, modulo diffeomorphism, on $n\choose2$ functions of $n$ variables. (For $n=2$, there is a slight difference in the numbers because $O(2)$ does not act with finite stabilizer on the generic $2$-jets, but it does act with finite stabilizer on the generic $k$-jets for $k>2$, so the answer to Question (1) is 'yes' in this case as well.)
As far as the answer to Question (2), for $S^{(k)}/\mathrm{Diff}_{k+1}$, it reduces to finding a description of the moduli space of orbits of the compact group $\mathrm{O}(n)$ on a finite dimensional vector space, which is a question in the invariant theory of representations of compact groups. (Generally, such a moduli space is embedded into some $\mathbb{R}^m$ as a closed semi-algebraic set by a set of $m$ generators of the ring of $\mathrm{O}(n)$-invariant polynomials on the vector space, but, for all but the simplest cases, finding such a set of generators and describing their relations is very difficult.)
Finally, answers for Question (3) are generally available. They depend on the theory of involutive systems of PDE, which I'm most familiar with in the form developed by Élie Cartan. For example, for hyper-Kähler metrics in dimension $4$ (which, interestingly, Cartan treated in 1926, long before the term 'hyper-Kähler' was invented), one finds that the $k$-jets of such metrics can be written in 'normal coordinates' in a manner analogous to that above for general metrics in such a way that they are the orbit space of $\mathrm{SU}(2)$ acting on a finite dimensional vector space $\mathcal{H}^k$. These vector spaces have the property that
$$
\dim \mathcal{H}^k - \dim \mathcal{H}^{k-1} = (k+3)(k-1)
$$
for $k\ge 2$, and $\mathcal{H}^{0}= \mathcal{H}^{1}= (0)$. The generic $\mathrm{SU}(2)$-stabilizer of a $2$-jet is finite, so the dimension of the $k$-jet orbit space is $\dim\mathcal{H}^k - 3$ for $k\ge2$. In particular, the analogous Poincaré series is indeed rational (which turns out to be a quite general phenomenon).
Remark: For more details on how this works, one might consult these notes of mine on exterior differential systems and their applications to problems of this kind.
As a particular example (just to keep the discussion simple), consider the case of describing the space of germs of the $G$-structures for $G\subset\mathrm{O(n)}$ that have vanishing intrinsic torsion and whose curvature tensor $R$ is required to satisfy some identity of the form $f(R) = 0$ (which may be trivial).
Some examples are (i) Riemannian metrics ($G=\mathrm{O}(n)$; no relation on the curvature), (ii) Ricci-flat metrics ($G=\mathrm{O}(n)$; $\mathrm{Ric}(R)=0$), (iii) Conformally flat metrics in dimensions $n\ge 4$ ($G=\mathrm{O}(n)$; $\mathrm{Weyl}(R)=0$), (iii) Self-dual Einstein metrics in dimension $n=4$ ($G=\mathrm{O}(4)$; $\mathrm{Weyl}_-(R)=\mathrm{Ric}^0(R)= 0$), (v) Kähler metrics ($G=\mathrm{U}(\tfrac12n)\subset\mathrm{O}(n)$; no relation on the curvature), (vi) torsion-free $\mathrm{G}_2$-structures $n=7$ ($G=\mathrm{G}_2\subset\mathrm{O}(7)$; no relation on the curvature), etc.
The possible curvature tensors for such a $G$-structure take values in a subset $A\subset \mathrm{Sym}^2\bigl(\Lambda^2(\mathbb{R}^n)\bigr)$ that is frequently a subspace invariant under the natural action of $G$ on $\mathrm{Sym}^2\bigl(\Lambda^2(\mathbb{R}^n)\bigr)$. In such cases, $A$ will often have the property that it is an involutive tableau of second order when regarded as a subspace of $\Lambda^2(\mathbb{R}^n)\otimes \Lambda^2(\mathbb{R}^n)$, with Cartan characters $s_2,\ldots,s_n$. Assuming that the PDE represented by $f(R)=0$ is, itself, involutive (again, true in all the cases listed above), one finds that the space of diffeomorphism classes of $k$-jets of $G$-structures satisfying $f(R)=0$ can be identified with the quotient (by the natural action of $G$) of the vector space
$$
A\oplus A^{(1)}\oplus \cdots\oplus A^{(k-2)}
$$
where $A^{(i)}$ is the $i^{th}$-prolongation of $A$ (as a second order tableau). Consequently, the Poincaré series of the graded vector space
$$
\mathcal{A} = A^{(0)}\oplus A^{(1)}\oplus \cdots\oplus A^{(k)}\cdots
$$
(where $A = A^{(0)}$ and where $A^{(i)}$ is assigned degree $i{+}2$) is
$$
P(t) = t^2\left(\frac{s_2}{(1{-}t)^2} + \cdots
+\frac{s_n}{(1{-}t)^n}\right).
$$
So, for example, for the general metric in dimension $n$, we find
$$
s_k = \tfrac12\,n(k{-}1)(n{-}k+1)\qquad 2\le k\le n
$$
while, for Ricci-flat metrics in dimension $n$, we find
$$
s_k = \tfrac12\,n(k{-}1)(n{-}k+1)\ \ (2\le k\le n{-}2);
\quad s_{n-1} = n(n{-}3);\quad s_n = 0.
$$
For the case of hyperKähler metrics in dimension $n=4m$, one finds that
$$
s_k = \tfrac12\,(k{-}1)(2m{+}2{-}k)(2m{+}3{-}k)\ \ (2\le k\le 2m{+}1)\ \text{and}
\ s_k = 0\ (2m{+}2\le k).
$$
Finally, for the case (vi) above (torsion-free $\mathrm{G}_2$-structures in dimension $7$), the system is involutive with
$$
(s_2,\ldots,s_7) = (14,21,21,15,6,0).
$$
