What does $\sum_{i_1=1}^n\sum_{i_2=1}^{i_1}\cdots \sum_{i_k=1}^{i_{k-1}} i_k^2$ equal? So, the figurate numbers $P_{k+1}(n) = \displaystyle \sum_{i_1=1}^n\sum_{i_2=1}^{i_1}\cdots \sum_{i_k=1}^{i_{k-1}}i_k = \binom{n+k}{k+1}$, which is quite neat.
I'm wondering what  the "square figurate numbers", $_{k+1}S(n)$, are equal to. This Wikipedia article gives the answer for $_2S(n)$, which is just $\displaystyle \frac{n(n+1)(2n+1)}6$. I'm looking for the general case, however.
$$_{k+1}S(n) = \sum_{i_1=1}^n\sum_{i_2=1}^{i_1}\cdots \sum_{i_k=1}^{i_{k-1}} i_k^2$$
Given that figurate number is used in a few different ways, see this Wikipedia article, under the section titled Triangular numbers and their analogs in higher dimensions, for the sense of the term that I'm using here.
 A: More generally $$\sum_{i_1=1}^n\sum_{i_2=1}^{i_1}\dots\sum_{i_k=1}^{i_{k-1}}f(i_k)=\sum_{i=1}^{n}c_{n,k,i}f(i),\qquad c_{n,k,i}=\binom{n-i+k-1}{k-1}$$ (indeed, $c_{n,k,i}$ counts sequences $(i_1,\dots,i_{k-1})$ such that $n\geqslant i_1\geqslant\dots\geqslant i_{k-1}\geqslant i$).
In terms of generating functions, if the above sum is $S(n,k,f)$, then $$\sum_{n=1}^\infty f(n)z^n=g(z)\implies\sum_{n=1}^\infty S(n,k,f)z^n=\frac{g(z)}{(1-z)^k}.$$
In our case $f(n)=n^2$ gives $g(z)=\dfrac{z(1+z)}{(1-z)^3}$ and the answer $\dbinom{n+k}{k+2}+\dbinom{n+k+1}{k+2}$.
A: First, a simple observation:
$$ n^2 = \underbrace{n+n+\ \dots \ + n}_{n \  \text{times}}$$
This observation allows us to switch over to a sigma function.
$$n^2 = n+ \sum_{i=1}^{n-1} (n-i) + i  = n + 2\binom n2$$
From here, evaluating $_{k+1}S(n)$ follows simply from the evaluation of $P_k(n)$.
$$\begin{align} _{k+1}S(n) &= \sum_{i_1=1}^n \sum_{i_2=1}^{i_1}\cdots \sum_{i_k=1}^{i_{k-1}} i_k^2 \\[2ex] &= \sum_{i_1=1}^n \sum_{i_2=1}^{i_1}\cdots \sum_{i_k=1}^{i_{k-1}}  i_k + 2\binom{i_k}2 
 \\[2ex] &=  \binom{n+k}{k+1} +  2\binom{n+k+1}{k+2}\end{align} $$
