Equivalent of Pauli matrices in 4 dimensions I would like to decompose the following 4x4 matrix:
$$ \mathrm{H} = 
    \begin{pmatrix}
    a & b & b & 0 \\
    b & 0 & 0 & b \\
    b & 0 & 0 & b \\
    0 & b & b & (-a+c)\\
    \end{pmatrix}
$$
in such a way that to computing the exponential of this matrix 
would have equivalent representation as the generalised Euler's formula
$$e^{ia(\hat{n}\cdot\vec{\sigma})} = \Bbb{1}\operatorname{cos}(a)+i(\hat{n}\cdot\vec{\sigma})\operatorname{sin}(x)\tag{1}\label{eq1}$$
with 
$$ \mathrm {M} = a(\hat{n}\cdot\vec{\sigma})  $$
M being the initial matrix.
Where $\vec{\sigma}$ is the so called Pauli vector containing the Pauli matrices as elements, and $\hat{n}$ is the normalised vector with coefficients constituting the decomposition of any 2x2 matrix regarding the Pauli matrices.
1 in the above represents the 2x2 dim unit matrix.
Is  there a an analogue to spin matrices in 4x4 dim, which can serve as the basis for this decomposition? 
 A: I'm a bit unclear as to what you actually want but let me point out a couple of things of what you have — I suspect you don't really mean to have c, since, unlike the Pauli vector, your matrix is not traceless: it has trace c.  I'll ignore it below, for simplicity. It doesn't modify the reduction of dimension discussed.
Your residual matrix 
$$ \mathrm{H} = 
    \begin{pmatrix}
    a & b & b & 0 \\
    b & 0 & 0 & b \\
    b & 0 & 0 & b \\
    0 & b & b & -a\\
    \end{pmatrix}
$$
is redundant, since it has the null eigenvector $(0,1,-1,0)$, so you may define 
an orthogonal transformation decoupling it,
$$  R = 
    \begin{pmatrix}
    1 & 0 & 0 & 0 \\
    0 & 1/\sqrt{2} &  1/\sqrt{2} & 0 \\
    0 &  -1/\sqrt{2} &  1/\sqrt{2} & 0 \\
    0 & b & b &  1\\
    \end{pmatrix}
$$
i.e. $R R^T=1\!\!1$,
$$
R H R^T=   \begin{pmatrix}
    a & \sqrt{2} b & 0 & 0 \\
    \sqrt{2} b & 0 & 0 & \sqrt{2}b \\
    0 & 0 & 0 & 0 \\
    0 & \sqrt{2}b & 0 & -a\\
    \end{pmatrix}
$$
so the 3rd component of your vector space is dross, and your your matrix is really just
$$
 \begin{pmatrix}
    a & \sqrt{2} b &  0 \\
     \sqrt{2}b & 0  & b\sqrt{2} \\
    0 & b\sqrt{2} &  -a\\
    \end{pmatrix}  =  a  \begin{pmatrix}
    1 &  0  & 0 \\ 0& 0 &0 \\
    0 & 0 &  -1\\
    \end{pmatrix} + b \sqrt{2} \begin{pmatrix}
    0 &  1 &  0 \\
     1 & 0  & 1 \\
    0 & 1 &  0\\
    \end{pmatrix} \equiv aA + bB.
$$
Note that, up to a normalization, A and B have the same eigenvalues, so they are orthogonal transforms of each other, and in that sense they are reminiscent of the two real, symmetric ones among the  Pauli matrices.
Since I am unclear about your problem (Your matrix is symmetric and so not an antisymmetric/anti-Hermitian generator of a rotation in 3D; it appears your invoking rotations might be misguided), I could suggest simply diagonalizing your matrix by an orthogonal transformation O, (its eigenvalues are 0 and $\pm \sqrt{a^2+2b^2}$), and simply exponentiating  the diagonal matrix. The orthogonal transform of the exponential of the diagonal matrix will be the exponential of the orthogonal transform of your matrix, and conversely. $O e^D ~ O^T= e^{ODO^T}.$
