I found a formula for counting 2-stars in a graph. Let $A$ be the adjacency matrix of a graph $\mathcal{G}$, then according to Snijders the number of 2-stars in is:
$$ \sum_{1 \leq i < j \leq g} \sum_{k \neq i,j} A_{i,k} \times A_{j,k} $$
where $g$ is the number of nodes in $\mathcal{G}$. This idea is easily extended to a three-stars:
$$ \sum_{1 \leq i < j < k \leq g} \sum_{l \neq i,j,k} A_{i,l} \times A_{j,l} \times A_{k,l} $$
and so on. Now continuing this approach should provide correct results, however it is computationally rather tedious, hence I wanted to ask if there is a more efficient way to calculate those numbers also for higher orders such as four-stars etc. Here is a simple python implementation:
def num_two_stars(A):
n = A.shape[0]
s = 0
for i in range(n):
for j in range(i+1, n):
for k in range(n):
if k != i and k!=j:
s += A[i,k] * A[j,k]
return s
def num_three_stars(A):
n = A.shape[0]
t = 0
for i in range(n):
for j in range(i+1, n):
for k in range(j+1, n):
for l in range(n):
if l != i and l!=j and l != k:
t += A[i,l] * A[j,l] * A[k,l]
return t
This is just a brute force algorithm $n^K$. You might find this paper useful. Also, you might possibly generate more interest in this on CS theory stack exchange (?)