Solution 1:

Maybe not exactly what you are asking, but notice that your matrix is lower triangular, and can be written as $M=I+N$, with $I$ the identity and \begin{equation} N=\begin{pmatrix} 0&0&0&0\\ 1&0&0&0\\1&1&0&0\\1&1&1&0 \end{pmatrix}.\end{equation}

N is nihilpotent, $N^4=0$ and (obviously) commutes with $I$, hence \begin{equation} (I+N)^n=\sum_{i=0}^n \binom{n}{i} N^i=\sum_{i=0}^3 \binom{n}{i} N^i, \end{equation} hence you only need $N$, $N^2$ and $N^3$ to compute any power of $M$.

Solution 2:

These are diagonals of Pascal's Triangle.

Solution 3:

We can simplify GFR's approach and go further: let $X$ be the matrix

$$ \left( \begin{matrix} 0&0&0&0\\1&0&0&0\\0&1&0&0\\0&0&1&0 \end{matrix} \right) $$

Then

$$M = I + X + X^2 + X^3 = \sum_{k=0}^{\infty} X^k = (I - X)^{-1}$$

Consequently,

$$M^n = (I-X)^{-n} = \sum_{k=0}^{\infty} \binom{-n}{k} (-X)^k = \sum_{k=0}^{\infty} \binom{n+k-1}{k} X^k = \sum_{k=0}^{3} \binom{n+k-1}{k} X^k $$

where I've used the generalized Binomial theorem and the formula for binomial coefficients to negative integers.