Midpoint iterative ellipse A set of ordered points $(x_i, y_i)$ on the plane are connected end to end, and each iteration takes the midpoint and connects to form a new polygon.
$$
(x_{k,i}, y_{k,i})=\left(
\frac{x_{k-1,i}+x_{k-1,i}}{2}, 
\frac{y_{k-1,i+1}+y_{k-1,i+1}}{2}
\right)
$$
The last point needs to be connected with the first point to calculate the midpoint:
$$
(x_{k,n}, y_{k,n})=\left(
\frac{x_{k-1,n}+x_{k-1,n}}{2}, 
\frac{y_{k-1,1}+y_{k-1,1}}{2}
\right)
$$
After enough iterations, the graph approximates an ellipse.

How to find the semi-major and semi-minor axes $a, b$ and the rotate angle $\alpha$ of this ellipse?
Is it possible to calculate the shape parameters of the final ellipse directly from the initial points without iteration?

Update 1
A similar transformation is required to zoom in when drawing, to prevent the precision from being reduced to zero, the code is as follows.
next[{xs_, ys_}] := Block[
    {x, y},
    x = ListConvolve[{1 / 2, 1 / 2}, xs, -1];
    y = ListConvolve[{1 / 2, 1 / 2}, ys, -1];
    (* Move to origin, won't change the shape *)
    x = x - Mean[x];
    y = y - Mean[y];
    (* Similarity transformation prevents exponential shrinking *)
    {x, y} / Max[Max[x] - Min[x], Max[y] - Min[y]]
];
drawPoints[this_] := Graphics[
    {
        PointSize[0.02], Blue, Point /@ this,
        Black, Line@Join[this, {First@this}]
    },
    PlotRange -> 0.6,
    ImageSize -> {300, 300}
];
drawAnimation[points_, nests_] := Block[
    {seq, w = 2, h = 1},
    seq = Transpose /@ NestList[next, {RandomReal[{-w, w}, points], RandomReal[{-h, h}, points]}, nests];
    drawPoints /@ Rest[seq] // ListAnimate
];

drawAnimation[25, 200]

 A: In certain sense, there is no final ellipse. The apparent ellipse will contract for a factor $\cos\frac{\pi}{n}$ in each iteration. At the end, all vertices will converge to a single point, the original vertex centroid of polygon.
Choose a coordinate system where the original vertex centroid is origin. Let $u_0 \in \mathbb{C}^n$ be an $n \times 1$ column vector with entries $z_k = x_k + y_k i$ at $k^{th}$ row. After $t$ iteration, we will
use $(z_k)_t$ and $u_t$ to denote the location of vertices and corresponding "u" vector.
In each iteration, the "$u$" at step $t$ and $t+1$ are related by a matrix equation
${}^{\color{blue}{[1]}}$
$$u_{t+1} = \Lambda u_t\quad\text{ where }\quad
\Lambda_{jk} = \begin{cases} \frac12, & k - j \equiv 0, 1 \pmod n\\0,& \text{ otherwise }\end{cases}.$$
Start from a polygon with vertex centroid at origin and repeat apply $\Lambda$ to "$u$". After enough number of iterations, $u_t$ will be dominated by eigenvectors of $\Lambda$ which are not orthogonal to original $u_0$ whose eigenvalues are largest in magnitude.
The largest eigenvalue of $\Lambda$ is $1$. The corresponding eigenvector $v_1 \propto (1,1,\cdots)^T$. Since we start with vertex centroid at origin, the original $u_0$ is orthogonal to $v_1$ and $v_1$ will not affect the asymptotic behavior of the polygon.
The next two largest eigenvalues in magnitude are degenerate.
They and corresponding eigenvectors have the form:
$$\lambda_{\pm} = \frac12(1 + \omega^{\pm 1})
\quad\longleftrightarrow\quad
v_{\pm} = (1, \omega^{\pm 1}, \omega^{\pm 2}, \cdots, \omega^{\pm(n-1)})^T
$$
where $\omega = e^{i\frac{2\pi}{n}}$.
If the vertices of the polygon are chosen randomly, then with probability $1$, the original $u_0$ will have non-zero projection on both $v_{\pm}$. So for large $t$, we have
$$u_t \sim A_{+} \lambda_{+}^t v_{+} + A_{-} \lambda_{-}^t v_{-}$$
for some coefficients $A_{\pm}$.
Since $|\lambda_{+}| = |\lambda_{-}| = \cos\frac{\pi}{n}$, you will find the polygon
contract for a factor $\cos\frac{\pi}{n}$ at each iteration.
With probability $1$ again, both $A_{\pm}$ will be non-zero. If you work out where the $z_k$ are positioned, you will find they lie on an ellipse. This explains why the polygon approaches an ellipse for large $t$.
To get the semi-major/semi-minor axes of the "ellipse" at large $t$, you need values of $A_{\pm}$. They are simply the "projection" of $u_0$ onto $v_{\pm}$. ie.
$$A_{\pm} = \frac1n \sum_{k=0}^{n-1} \omega^{\mp k} z_k $$
I will leave the actual computation of semi-major/semi-minor axes to you.
Note

*

*$\color{blue}{[1]}$  - vertices are labelled such that $(z_k)_{t+1} = \frac12(z_k + z_{k+1})_t$. ie. the $k^{th}$ vertex at step $t+1$ is the mid-point of the segment joining $k^{th}$ and $(k+1)^{th}$ vertices at step $t$.

A: For centroid origin and after sufficient iterations, the curve can be fitted approximately by the conic
\begin{align}
  bx^2-2hxy+ay^2 &= 2(ab-h^2)\\
  a &= \mathrm{Var}(X) \\
  b &= \mathrm{Var}(Y) \\
  h &= \mathrm{Cov}(X,Y)
\end{align}
Note that $$\mathrm{Var}(X) \mathrm{Var}(Y) \ge \mathrm{Cov}(X,Y) \implies ab-h^2 \ge 0$$
The above result will be invariant for iterating Steiner in-ellipse:
$$\left| z-2\sqrt{uv}\cos \frac{k\pi}{n} \right|+
\left| z+2\sqrt{uv}\cos \frac{k\pi}{n} \right|=
2(|u|+|v|) \left| \cos \frac{k\pi}{n} \right|$$
of $\{ \frac{n}{k} \}$ star-polygon with vertices $z_j=u\, \omega_n^j+v\, \omega_n^{(n-1)j}$ where $\omega_n=e^{\frac{2\pi i}{n}}$ and $u,v\in \mathbb{C}$.
