Why does $\cos(x) + \cos(y) - \cos(x + y) = 0$ look like an ellipse?

It looks like an ellipse because the isocurves of any smooth surface look like an ellipse in the vicinity of an extremum !

Indeed, by Taylor's development in 2D,

$$f(x,y)=f(a,b)+\frac{\partial f}{\partial x}(x-a)+\frac{\partial f}{\partial y}(y-b)\\ +\frac12\frac{\partial^2 f}{\partial x^2}(x-a)^2+\frac{\partial^2 f}{\partial x\partial y}(x-a)(y-b)+\frac12\frac{\partial^2 f}{\partial y^2}(y-b)^2\cdots$$

This development shows that in the vicinity of an ordinary point, a smooth surface usually behaves like a plane because the linear terms dominate, and isocurves are approximately straight lines.

But when the first derivatives are zero, the next order terms, the quadratic ones, enter into play and the behavior becomes that of a quadric. The isocurve $f(x,y)=F$ is approximately of the form

$$A(x-a)^2+2B(x-a)(y-b)+C(y-b)^2=D,$$ an ellipse centered at $(a,b)$. (Provided that the quadratic form is positive-definite, i.e. not a saddle point.)

The closer you get to the extremum, the more you get an exact ellipse. The effect is more pronounced and the curves more symmetric when the third order derivatives are small.

In the case of the given trigonometric function, translating the coordinates to $(\pi,-\pi)$, to bring a maximum at the origin, we get

$$z=\cos(x)+\cos(y)+\cos(x+y).$$

Then the Taylor development up to fourth order yields

$$z\approx 3-x^2-xy-y^2+\frac{x^4}{12}+\frac{x^3y}6+\frac{x^2y^2}4+\frac{xy^3}6+\frac{y^4}{12}.$$

The plot below shows you the isocurve $z=0$ computed with the quadric (brownish) and quartic (pinkish) approximations; the latter is indistinguishable from the true curve, a quasi-ellipse.

Using the sum-to-product formula for cosine

$$ \cos s + \cos t = 2\cos \tfrac{s+t}{2}\cos \tfrac{s-t}{2} = \cos (s+t)$$

this is a great time to do that 45 degree rotation $u = \tfrac{s+t}{2}:, v=\tfrac{s-t}{2}, s+t = 2u$

$$ 2\cos u\cos v = \cos 2u \hspace{0.25in}\text{or}\hspace{0.25in} \bbox[5px,border:2px solid #F5A029]{2\cos v = \frac{\cos 2u}{\cos u}}$$

Wolfram Alpha doens't offer much simplification beyond this.

import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline

x = 2*np.pi*np.arange(-1,1,0.01)
y = 2*np.pi*np.arange(-1,1,0.01)
z = np.cos(x[...,None]) + np.cos(y[None,...]) - np.cos(x[...,None]+y[None,...])
plt.contour(x,y,z)

They are definitely not ellipses:

Alternatives plots from Wolfram Alpha. Somehow the level sets must interpolate between the "ellipses" and the "bow-ties".

Amoebas and Coamoebas

Using Demoivre identity $\cos \theta + i \sin \theta = e^{i\theta}$ we can instead consider the equation in complex numbers:

$$ e^{ix} + e^{iy} - e^{i(x+y)} = z + w - zw = 1 - (z-1)(w-1) = 0$$

The set of points $(z,w) \in \mathbb{C}^2$ satisfying this equation is not an algebraic variety. The log of the norms are what is known as an amoeba.

$$ \log: (z,w) \mapsto (\log |z|, \log|w|) \hspace{0.25in}\text{for}\hspace{0.25in} (z,w) \in \{ (z,w):z + w - zw\} $$

The angles of the solutions to this equation are known as the coamoeba

$$ \arg: (|z|e^{i\theta},|w|e^{i\phi}) \mapsto (\theta, \phi) \hspace{0.25in}\text{for}\hspace{0.25in} (z,w) \in \{ (z,w):z + w - zw\} $$

See What is an Amoeba in the Notices of the American Mathematical Society. Or the Master's thesis of Masahito Yamazaki Brane Tilings and their Applications

Your equation is not an ellipse, or even a family of ellipses. Ellipses are graphs given by a quadratic relation in $$Ax^2+Bxy+Cy^2+Dx+Ex+F=0.$$ Your equation, using Taylor Expansion (for $k\to\infty$), is $$\sum _{n=0}^k\frac{\left(-1\right)^n\left(x^{2n}+y^{2n}\right)}{\left(2n\right)!}=\sum _{n=0}^k\frac{\left(-1\right)^n\left(x+y\right)^{2n}}{\left(2n\right)!},$$ which is only quadratic for $k=1,$ and incidentally the graph for $k=1$ is not an ellipse but the hyperbola $xy+1=0.$ As a related note, a near duplicate of your graph using an ellipse is $$\left(x-\pi \right)^2+\left(x-\pi \right)\left(y-\pi \right)+\left(y-\pi \right)^2-\sqrt{19}=0.$$ But by inspection we see the curvature of the graphs are not equal.

In fact, here is an animation showing the polynomial relations for $k=0$ to $k=20.$

The "ellipses" occur along the lines $y=\pm x.$

Also, as conjectured my @MathGemini in the above comments, your graph could be an agglomeration of cones in $(x,y,z),$ where $$z=\cos x+\cos y-\cos(x+y).$$ This is not the case, however because the surface you get is one of periodic maxima and minima, which is the expected behavior given by the equation.

(This is the surface and its intersection with the $xy$ plane viewed from the positive $z-$axis.)

(This is the surface and its intersection with the $xy$ plane viewed from $x>0,y>0,z>0.$)

(This is the surface and its intersection with the $xy$ plane viewed from below the $xy-$plane.)

They are not ellipses.

To prove this, consider the closed curve lying above and directly to the right of the origin. Since the given equation is symmetric between $x$ and $y$, this curve is symmetric about the line $y=x$. Thus, if it is an ellipse, its minor axis must line on the line $y=x$. Now, it is easy to check that the points $$ (\pi,\pi) \pm \alpha (1,1) \qquad\text{and}\qquad (\pi,\pi) \pm \beta (-1,1), $$ lie on this curve, where $\alpha = \cos^{-1}\!\!\left(\dfrac{-1+\sqrt{3}}{2}\right)$ and $\beta = \dfrac{2\pi}{3}$, so these would have to be the endpoints of the major and minor axes. However, in this case, the point $$ (x,y) \;=\; (\pi,\pi) + \frac{\alpha}{\sqrt{2}}(1,1) + \frac{\beta}{\sqrt{2}}(-1,1) $$ would also have to lie on the ellipse, and it doesn't. In particular, $$ \cos(x) + \cos(y) -\cos(x+y) \;\approx\; -1.906 \;\ne\; 0. $$ for these values of $x$ and $y$.

For the record, the substitution $t=\cos x, s=\cos y$ leads to $2st(s+t) + 1 = 2(s^2+t^2+st)$ with $s, t$ ranging over $[-1, 1]$. We now need some algebraic geometer to recognize this equation...