Why not use the third derivative for numerical optimization? If Hessians are so good for optimization (see e.g. Newton's method), why stop there? Let's use the third, fourth, fifth, and sixth derivatives? Why not?
 A: Even calculating Hessians is quite a bit of work:
$$H = \begin{bmatrix}
  \dfrac{\partial^2 f}{\partial x_1^2} & \dfrac{\partial^2 f}{\partial x_1\,\partial x_2} & \cdots & \dfrac{\partial^2 f}{\partial x_1\,\partial x_n} \\[2.2ex]
  \dfrac{\partial^2 f}{\partial x_2\,\partial x_1} & \dfrac{\partial^2 f}{\partial x_2^2} & \cdots & \dfrac{\partial^2 f}{\partial x_2\,\partial x_n} \\[2.2ex]
  \vdots & \vdots & \ddots & \vdots \\[2.2ex]
  \dfrac{\partial^2 f}{\partial x_n\,\partial x_1} & \dfrac{\partial^2 f}{\partial x_n\,\partial x_2} & \cdots & \dfrac{\partial^2 f}{\partial x_n^2}
\end{bmatrix}.$$
Now see how the third derivative looks like:
$$\partial H/\partial x=\begin{bmatrix}
  \dfrac{\partial H}{\partial x_1}\\
  \dfrac{\partial H}{\partial x_2}\\
  \vdots\\
  \dfrac{\partial H}{\partial x_n}
\end{bmatrix}$$
This is a three dimensional matrix. Here's how its elements look like:
$$(\partial H/\partial x)_{ijk}=\dfrac{\partial^3 f}{\partial x_i\partial x_j\partial x_k}$$
Sixth's derivative will be six dimensional matrix:
$$\dfrac{\partial^6 f}{\partial x_i\partial x_j\partial x_k\partial x_l\partial x_m\partial x_n}$$
Usually, the trade-off is not favorable for going after higher than Hessian. I mean the trade-off between potential gain in speed through using higher order approximations vs. the noise amplification. You always have noise in inputs because we're talking about statistical applications. This noise will be amplified by the derivatives. 
If you play golf then the analogy in optimization is to first swing trying to get to the green, not worry to much about a hole. Once, on the green, we'll putt aiming a hole. 
A: Typically, when you analyze the effectiveness of such algorithms, you'll find results such as one step of a fourth order algorithm having roughly the same effectiveness as two steps of a second order algorithm.
So the choice of which algorithm to use is relatively simple: if one step of the fourth order algorithm takes twice as much work or more than one step of the second order algorithm, you should use the latter instead.
That is the typical situation for these sorts of methods: the classical algorithm has the optimal work-to-effectiveness ratio for general problems. While there are occasional problems where a higher order approach is unusually easy to compute and can outperform the classical variant, they are relatively uncommon.
A: I am interpreting the question as being "Why does Newton's method only use first and second derivatives, not third or higher derivatives?"
Actually, in many cases, going to the third derivative does help; I've done it with custom stuff before.  However, in general, going to higher derivatives adds computational complexity - you have to find and calculate all those derivatives, and for multivariate problems, there are a lot more third derivatives than there are first derivatives! - that far outweighs the savings in step count you get, if any.  For example, if I have a 3-dimensional problem, I have 3 first derivatives, 6 second derivatives, and 10 third derivatives, so going to a third-order version more than doubles the number of evaluations I have to do (from 9 to 19), not to mention increased complexity of calculating the step direction / size once I've done those evaluations, but will almost certainly not cut the number of steps I have to take in half.
Now, in the general case with $k$ variables, the collection of $n^{th}$ partial derivatives will number${k+n-1} \choose {k-1}$, so for a problem with five variables, the total number of third, fourth, and fifth partial derivatives will equal 231, a more than 10-fold increase over the number of first and second partial derivatives (20).  You would have to have a problem that is very, very close to a fifth-order polynomial in the variables to see a large enough reduction in iteration counts to make up for that extra computational burden.
A: You can think of the order of derivatives as the order of a polynomial approximation to the function. Most optimization routines rely on convexity. A quadratic polynomial will be convex/concave everywhere whereas a 3rd order or higher polynomial will not be convex everywhere. Most optimization routines rely on successive approximations of convex functions with quadratics for this reason. A quadratic approximation that is convex requires a positive definiteness condition to be imposed in order for the quadratic to be convex. 
A: I don't really see what the statistical aspect of this question is, so I'll answer the optimization part.
There are 2 parts to convergence: iteration cost & iteration count
Pretty much every answer here is focusing on just the iteration cost and ignoring the iteration count. But both of them matter. An method that iterates in 1 nanosecond but that takes $10^{20}$ iterations to converge won't do you any good. And a method that blows up won't help either, no matter how cheap its iteration cost.
Let's figure out what's going on.
So: Why not use >2nd-order derivatives?
Partly because (and this is true for 2nd-order too, but more on that in a bit):
Higher-order methods generally only converge faster when near the optimum.
On the other hand, they blow up more easily when they are farther from the optimum!
(Of course, this isn't always true; e.g. a quadratic will converge in 1 step with Newton's method. But for arbitrary functions in the real world that don't have nice properties, this is generally true.)
This means that when you are farther away from the optimum, you generally want a low-order (read: first-order) method. Only when you are close do you want to increase the order of the method.
So why stop at 2nd order when you are near the root?
Because "quadratic" convergence behavior really is "good enough"!
To see why, you first have to understand what "quadratic convergence" means.
Mathematically, quadratic convergence means that, if $\epsilon_k$ is your error at iteration $k$, then the following eventually holds true for some constant $c$:
$$\lvert\epsilon_{k+1}\rvert \leq c\ \lvert\epsilon_{k}\rvert^2$$
In plain English, this means that, once you are near the optimum (important!), every extra step doubles the number of digits of accuracy.
Why? It's easy to see with an example: for $c = 1$ and $\lvert\epsilon_1\rvert = 0.1$, you have $\lvert\epsilon_2\rvert \leq 0.01$, $\lvert\epsilon_3\rvert \leq 0.0001$, etc. which is ridiculously fast. (It's super-exponential!)
Why not stop at 1st order rather than 2nd-order?
Actually, people often do this when second-order derivatives become too expensive. But linear convergence can be very slow. e.g. if you got $\epsilon_k = 0.9999999$ then you'd need maybe 10,000,000 iterations with linear convergence to get $\lvert\epsilon\rvert < 0.5$, but only 23 iterations with quadratic convergence. So you can see why there's a drastic difference between linear and quadratic convergence. This is not true for 2nd and 3rd-order convergence, for example (see next paragraph).
At this point, if you know any computer science, you understand that with 2nd-order convergence, the problem is already solved. If you don't see why, here's why: there is nothing practical to gain from tripling the number of digits every iteration instead of doubling it—what's it going to buy you? After all, in a computer, even a double-precision number has 52 bits of precision, which is around 16 decimal digits.
Maybe it will decrease the number of steps you require from 16 to 3... which sounds great, until you realize it comes at the price of having to compute third derivatives at each iteration, which is where the curse of dimensionality hits you hard. For a $6$-dimensional problem, you just paid a factor of $6$ to gain a factor of $\approx 5$, which is dumb. And in the real world problems have at least hundreds of dimensions (or even thousands or even millions), not merely $6$! So you gain a factor of maybe 20 by paying a factor of, say, 20,000... hardly a wise trade-off.
But again: remember the curse of dimensionality is half the story.
The other half is that you generally get worse behavior when you're far from the optimum, which generally adversely affects the number of iterations you have to do.
Conclusion
In a general setting, higher-order methods than 2 are a bad idea. Of course, if you can bring additional helpful assumptions to the table (e.g. perhaps your data does resemble a high-degree polynomial, or you have ways of bounding the location of the optimum, etc.), then maybe you can find that they are a good idea—but that will be a problem-specific decision, and not a general rule of thumb to live by.
A: Let me be the only one here defending 3rd order methods for SGD convergence, but definitely not in the entire space what would need $\approx dim^3/6$ coefficients, but e.g. in just a single direction, which needs only a single additional coefficient if already having 2nd order model in this direction.
Why single direction 3rd order model can be beneficial? For example because close to zero second derivative in this direction basically means two alternative scenarios: plateau or inflection point - only the former requires larger step size, and 3rd derivative allows to distinguish them.
I believe we will go toward hybrid multi-order methods: 2nd order method in a low dimensional subspace e.g. from PCA of recent gradients, what still allows for free 1st order simultaneous gradient descent toward part of gradient orthogonal to this subspace ... and additionally I would add e.g. 3rd order model for a single most relevant direction.
