Great question, but I'm not sure a good answer exists (at this point in time). Indeed, this estimator (alternatively called the REINFORCE estimator, the score function estimator [SFE], and the likelihood ratio estimator) is known to have very high variance, which is a major problem in RL, as well as in other problems (e.g., differentiating through discrete latent variable models).
I think there are three reasons, intuitively, why it has very high variance.
We use only values of $f$, never its derivative (which is usually unknown or non-existent). In other words, the estimator has no access to any information about how $f$ varies locally. Presumably, if we knew how our alteration was going to change $f$ infinitesimally, we could do better; i.e., exactly how our change in $\psi$ would affect $\theta$ (which we know probabilistically), and how perturbing $\theta$ would affect $f$. (We could then basically invoke the chain rule). Basically, we just don't have much information. Since we don't know how $f$ changes, we cannot build it into the estimator (and our estimator is of a gradient, i.e., something that is measuring how the function is supposed to be changing).
The SFE is extremely general, requiring absolutely nothing from $f$ except the need to be able to evaluate it. I think it is intuitive that it must pay a price for such generality. This is unlike, e.g., the reparametrization trick, which requires certain properties from $f$ and/or $\theta$, but may have lower variance (e.g., see ).
The SFE is unbiased. There are closely related methods that sacrifice this property for greatly reduced variance. I think it's intuitive that there is some bias-variance tradeoff for such estimators. Building off stochastic optimal control theory, it turns out it's possible to use control variate baselines to reduce the variance without introducing bias however (e.g., ).
The main alternative to the SFE is the reparameterization trick (closely related to pathwise derivative estimators, e.g. see [3,4]). However, if $f$ is truly a blackbox function, it cannot readily be applied (unlike the SFE). Furthermore, backprop through discrete variables requires "softening" the variable (e.g., the concrete  or Gumbel-Softmax  method), meaning some bias can be introduced to the estimator. Nevertheless, the resulting estimator is empirically known to be much lower variance and more stable (i.e., in applications).
I say empirically because it is known that theoretically the SFE can have lower variance in some cases (see , page 34). The same reference gives a condition (in the Gaussian case) under which the SFE is worse. In other words, it is theoretically possible that (in some cases) the SFE has lower variance than its alternatives, but in practice this usually seems to not be the case (again also see ).
 Variance reduction properties of the reparameterization trick, Xu et al, 2019
 Backpropagation through the void: Optimizing control variates for black-box gradient estimation, Grathwohl et al, 2017
 Pathwise Derivatives Beyond the Reparameterization Trick, Jankowiak & Obermeyer, 2018
 Implicit Reparameterization Gradients, Figurnov et al, 2018
 The concrete distribution: A continuous relaxation of discrete random variables, Maddison et al, 2016
 Categorical reparameterization with gumbel-softmax, Jang et al, 2016 (abstract)
 Uncertainty in Deep Learning (thesis), Gal, 2016
This is a rather active research area. See also: