The key thing here is that Wooldridge makes the assumption of random sampling.
Notice that since we have a random sample this means $(x_i, y_i) \perp (x_j, y_j)$ for $i \neq j$ which means that the single components of the pair are also independent, in particular $x_i \perp x_j$ (to see that just notice the joint is $p(x_i, y_i, x_j, y_j) = p(x_i, y_i)p (x_j, y_j)$ and marginalize over $y$).
This further implies $y_i \perp y_j |x_i, x_j$, since:
$$
p(y_i, y_j|x_i, x_j) = \frac{p(x_i, y_i, x_j, y_j)}{p(x_i, x_j)} = \frac{p(x_i, y_i)p (x_j, y_j)}{p(x_i)p(x_j)} = p(y_i|x_i)p(y_j|x_j)
$$
But $y_i|x_i$ is nothing more than the disturbance $u_i$ plus a constant. Hence $u_i \perp u_j |x$ when assuming a random sample and $E[u_i u_j | x] =E[u_i|x] E[u_j | x] = 0$
A harmonic mean is the reciprocal of the mean reciprocal of data, and is used to average rates. For example, electrical capacitance of a series of $n$ connected capacitors is $\frac{1}{n}^{th}$ of the harmonic mean, and electrical resistance of $n$ parallel connected resistors is also $\frac{1}{n}^{th}$ of its harmonic mean. In other words, for resistors in parallel, the harmonic mean resistance is the mean value of each individual resistance. Here, we use $t$ for time, but the following is true for any $x$. For continuous density functions one uses a variation of the second mean value theorem for integrals, which for support on $[\alpha\geq0,\beta]$, where $\alpha<\beta$ and where $1=\int_{\alpha }^{\beta } f(t) \, dt$, the harmonic mean residence time (H-MRT) is
\begin{equation}
\text{H-MRT} := \Bigg \langle \frac{1}{T} \Bigg \rangle ^{-1}=\Bigg \langle \frac{1}{T} \Bigg \rangle ^{-1}\Bigg[\int_{\alpha }^{\beta } f(t) \, dt\Bigg]^{-1}=\Bigg[\int_{\alpha }^{\beta } \frac{f(t)}{t} \, dt\Bigg]^{-1}.
\end{equation}
For example, the Pareto density, whose first moment is undefined (unphysical; negative) for $0<\alpha<1$, some measure other than the mean, e.g., the harmonic mean is needed when the right tail is very heavy. From the definition above, this is
\begin{equation}
\text{H-MRT}_{\text{Pareto}} := \Bigg[\int_{\beta }^{\infty } \frac{\alpha \beta ^{\alpha } t^{-\alpha -1}}{t} \, dt\Bigg]^{-1} = \beta \Big(1+\frac{1}{\alpha }\Big);\; \alpha,\beta>0
\end{equation}
For the case when the mean value of a Pareto distribution is undefined, three other statistical measures are defined; median, geometric mean, and the harmonic mean, all three of which are also defined when the first moment is also defined. Of those three, the harmonic mean is arguably most useful as argued in https://arxiv.org/pdf/1402.4061; "However, the harmonic mean statistic will be relatively insensitive to the value of $\alpha _{\min }$ and will tolerate an $\alpha _{\min }$ that is close to 0. That is another point in favor of using the harmonic mean instead of one of the other statistics."
For a left truncated exponential distribution, the harmonic mean is
\begin{equation}
\Bigg[ \int_{\alpha }^{\infty } \frac{\lambda e^{-\lambda(x-\alpha)}}{x} \, dx \Bigg]^{-1}=\frac{e^{-\alpha \lambda }}{\lambda \Gamma (0,\alpha \lambda )};\; \alpha >0,\lambda >0
\end{equation}
The left and right truncated exponential harmonic mean is
\begin{equation}
\Bigg[\int_{\alpha }^{\beta } \frac{\lambda e^{\lambda (\beta -x)}}{x \left(e^{\lambda (\beta -\alpha )}-1\right)} \, dx \Bigg]^{-1} =\frac{e^{-\beta \lambda } \left(e^{\lambda (\beta -\alpha )}-1\right)}{\lambda (\text{Ei}(-\beta \lambda )-\text{Ei}(-\alpha \lambda ))};\alpha <\beta ,\alpha \neq 0,
\end{equation}
where the exponential integral function $\text{Ei}(z)=-\int_{-z}^{\infty } \frac{e^{-t}}{t} \, dt$, where the principal value of the integral is taken.
For the beta distribution the harmonic mean is
\begin{equation}
\Bigg[ \int_0^1 \frac{x^{\alpha -1} (1-x)^{\beta -1}}{x B(\alpha ,\beta )} \, dx\Bigg]^{-1} =\frac{a-1}{a+b-1};\;\alpha>1,\beta>0
\end{equation}
Best Answer
You could point out that this is a weighted least squares regression with weights $1/x_i$.
To make the connection with the references, revert to a standard notation in which you seek to find $\beta$ that minimizes $$\sum \omega_i (y_i - \beta)^2.$$
This is a model with a single constant regressor $$X = \pmatrix{1\\1\\\vdots\\1}$$ and weights matrix $$W = \pmatrix{\omega_1 & 0 & \cdots & 0 \\ 0 & \omega_2 & \cdots & 0 \\\vdots & \vdots & \ddots & 0 \\ 0 & \cdots & 0 & \omega_n}.$$
I have renamed "$x_i$" as "$y_i$" (the "response") and the parameter to be estimated is $\beta$ instead of $z$. The weights are $\omega_i=1/x_i$. It is necessary that they all exceed $0$. The solution is
$$\hat\beta = (X^\prime W X)^{-1}X^\prime W y = \frac{\sum_i x_i\omega_i }{\sum_i \omega_i} = \frac{\sum_i x_i/x_i }{\sum_i 1/x_i} = \frac{n}{\sum 1/x_i},$$
QED.
Comments
The same analysis applies to any positive sets of weights, providing a generalization of the harmonic mean and a useful way to characterize it.
When, as in a controlled experiment, the $x_i$ are viewed as fixed (and not random), the machinery of weighted least squares provides confidence intervals and prediction intervals, etc. In other words, casting the problem into this setting automatically gives you a way to assess the precision of the harmonic mean.
Viewing the harmonic mean as the solution to a weighted problem provides insight into its nature and, especially, to its sensitivity to the data. It is now clear that the most important contributors are those with the smallest values of $x_i$--and their importance has been quantified by the weights matrix $W$.
Reference
Douglas C. Montgomery, Elizabeth A. Peck, and G. Geoffrey Vining, Introduction to Linear Regression Analysis. Fifth Edition. J. Wiley, 2012. Section 5.5.2.