Revised: From the clarifications it appears that the conical approximation of data points is only a rough one. However the apex is known, which leads to a revised suggestion of projecting points onto a sphere centered at the apex, which would lead to a rough circle of points on the sphere (these points would be exactly on a "small" circle of the sphere if the original data were actually on a cone).
Given the points on a sphere, one might proceed by fitting a small circle (great circles on a sphere being a special case) in a total least squares sense, i.e. minimizing the sum of squares of spherical distances from points to the fitted circle. This is a nonlinear problem, but it involves only three parameters (the radius and two coordinates for center of the circle). See Fujiki and Akaho, 2007 for exposition of both the 2D and higher dimensional analogs, and its references to the literature.
Your interest seems primarily to be identification of the axis of rotation, i.e. the pole that corresponds to the spherical circle. See here for some software tailored to this problem.
I think we should be careful about an approach that assumes points having a roughly uniform distribution over (around) the surface of a cone. If they were, it would make spotting a "V" shape in a slice of the data more likely, but perhaps there are more robust approaches (less dependent on an even sampling of points).
[Something else to keep in mind is that mathematically a cone proceeds in both directions of the axis through its apex. Therefore an ideal slice of the data through a plane containing the apex might reveal an "X" rather than a "V" shape.]
One possibility would be to assess for a candidate "apex" point $A$ of a cone, what axis gives the most tightly clustered slopes of lines passing through the sample points and the "apex". The axis is generated by a unit vector $U$ from $A$, and the slopes through sample point $P$ can be associated with a dot product of $U$ with $P - A$ divided by $||P - A||$.
In fact we might weight these slopes so that samples close to $A$ are less important than points far away (since the points close to $A$ are close to the cone regardless of what axis is chosen), and this suggests weighting by $||P - A||$, thereby cancelling the division used above.
This replaces the search for a "good slice" of the data with a search for a good approximation of the cone's apex.
The problem is of some importance in the field of computer vision and robotics. In this book chapter by Chernov and Ma, Least Squares Fitting of Quadratic Curves and Surfaces (2011), some improvements are sought using a "fast projection" method of Dave Eberly. The authors treat three non-degenerate quadratic surfaces in 3D, ellipsoids, hyperbolic paraboloids, and hyperboloids of one sheet, and mention without details that elliptic paraboloids and hyperboloids of two sheets can be treated similarly. A cone, however, is a degenerate quadric surface, intermediate between the hyperboloids of one and two sheets. I'd expect that the procedures described by Chernov and Ma would adapt readily to this specialization.
Best Answer
I think I found a solution.
The point $(x,y,z)$ is not only on the hyperbola, but also on the circle that results from intersecting the cone with the plane defined by $x$ and that is given by the equation $y^2+z^2=r^2$, where $r=\sqrt{y^2+z^2}$ is the radius of the circle, which has its center in $(x, 0, 0)$.
Now that I have the radius $r$ and length of the cone $x$, I can find the angle of the cone $\alpha=arctan(r/x)$, which was what I was looking for.
Tl;dr: $\alpha=arctan(\sqrt{y^2+z^2}/x)$.