I am filling in the details of the calculation for the Hessian, as we can derive a couple of relations that will save some work further along. You presumably found the critical points from
$$ yz \ = \ 2x \ \ , \ \ xz \ = \ 2y \ \ , \ \ xy \ = \ 2z \ \ \Rightarrow \ \ xyz^2 \ = \ 4xy \ = \ 8z \ \ ; $$
[multiplying the first two equations together, and also applying the third]
we may either take $ \ x \ = \ y \ = \ z \ = 0 \ $ , or, using this last equation with $ \ z \ \ne \ 0 \ $ ,
$$ xyz \ = \ 8 \ \ \Rightarrow \ \ xy \ = \ \frac{8}{z} \ = \ 2z \ \ \Rightarrow \ \ z^2 \ = \ 4 \ \ . $$
[applying the third critical point equation again]
From this, you found that $ \ z \ = \ \pm 2 \ $ , and then the coordinates you list.
The second partial derivatives you show produce the Hessian matrix (as you would have obtained)
$$ \left[ \begin{array}{ccc}-2&z&y\\z&-2&x\\y&x&-2\end{array} \right] \ \ , $$
for which the characteristic equation is
$$ - (2 + \lambda)^3 \ + \ 2 \ xyz \ - \ (x^2 + y^2 + z^2) \ [ -(2 + \lambda) \ ] \ = \ 0 \ \ . $$
The point of setting this up is that we can find first (as you did) that for the critical point at the origin, all of the eigenvalues are $ \ -2 \ $ (which also comes directly from the fact that the Hessian is diagonal in that case), but also that for all four of the other critical points,
$$ xyz \ = \ 8 \ \ , \ \ x^2 \ = \ y^2 \ = \ z^2 \ = \ 4 \ \ \Rightarrow \ \ - (2 + \lambda)^3 \ + \ 2 \cdot 8 \ - \ 3 \cdot 4 \ [ \ -(2 + \lambda) \ ] \ = \ 0 $$
$$ \Rightarrow \ \ \lambda^3 \ + \ 6 \ \lambda^2 \ - \ 32 \ = \ 0 \ \ \Rightarrow \ \ \lambda \ = \ 2 \ , \ -4 \ , \ -4 \ \ . $$
Since there are both positive and negative eigenvalues, all of the remaining critical points are saddle points [as Hayden remarks].
$$ \ $$
Another approach we can take is to look at "cross-sections" of the hypersurface described by $ \ f(x,y,z) \ $ near these critical points (in the absence of a convenient way to "visualize" the function). This turns out to be important because a simple assessment of values of the function can be deceptive. We have
$$ f(2, \ 2, \ 2) \ = \ f(-2, \ -2, \ 2) \ = \ f(-2, \ 2, \ -2) \ = \ f(2, \ -2, \ -2) $$
$$ = \ 8 \ - \ 3 \ (\pm 2)^2 \ = \ -4 \ \ . $$
For example, if we look at $ \ (2, \ 2, 2) \ $ and calculate values in the immediate vicinity of that points, we find (for instance)
$$ f(2.1, \ 2, \ 2) \ = \ f(1.9, \ 2, \ 2) \ = \ -4.01 \ \ ; $$
owing to the symmetry of the function, a similar small step in the $ \ y- \ $ or $ \ z-$ coordinate gives the same values. So one might reasonably form the impression that this point is a relative maximum, as the origin is.
If we take the "slice" $ \ x \ = \ y \ = \ z \ $ , however, our function becomes
$ \ f(x, \ x, \ x) \ = \ x^3 \ - \ 3x^2 \ $ , for which we find that $ (0, \ 0, \ 0) \ $ is a relative maximum (as we would expect) and $ \ (2, \ 2, \ 2) \ $ is a relative minimum (concavity is "upward" there). On the other hand, if we look at a "slice" perpendicular to that along $ \ x \ = \ z \ \ , \ \ y \ = \ 4 - x \ $ , the behavior of our function is $ \ f(x, \ 4 - x, \ x) \ = \ -x^3 \ + \ x^2 \ + \ 8x \ - \ 16 \ $ . Along this slice, $ \ (2, \ 2, \ 2) \ $ is a relative maximum ("downward" concavity). So this point is not a relative extremum of $ \ f(x,y,z) \ $ , but rather is a saddle point. Similiar analyses can be carried out for each of the other three critical points in question.