How about the following short "manual calculation"?

def weighted_avg_and_std(values, weights):
    """
    Return the weighted average and standard deviation.

    They weights are in effect first normalized so that they 
    sum to 1 (and so they must not all be 0).

    values, weights -- NumPy ndarrays with the same shape.
    """
    average = numpy.average(values, weights=weights)
    # Fast and numerically precise:
    variance = numpy.average((values-average)**2, weights=weights)
    return (average, math.sqrt(variance))

    

Interestingly, there is no single equation for a weighted standard error. Multiple versions have been proposed in the literature though. See for example:

Donald F. Gatz and Luther Smith (1995). "The Standard Error of a Weighted Mean Concentration - I: Bootstrapping Vs Other Methods". In: Atmospheric Environment 29.11, pp. 1185-1193

I implemented some of these method in R to use as an (unexported) function in the adjustedCurves R-package I developed, here is the code:

weighted.se <- function(x, w, se_method, na.rm=FALSE) {

  if (na.rm) {
    miss_ind <- !is.na(x)
    w <- w[miss_ind]
    x <- x[miss_ind]
  }

  n <- length(x)
  mean_Xw <- stats::weighted.mean(x=x, w=w, na.rm=na.rm)

  ## Miller (1977)
  if (se_method=="miller") {
    se <- 1/n * (1/sum(w)) * sum(w * (x - mean_Xw)^2)
  ## Galloway et al. (1984)
  } else if (se_method=="galloway") {
    se <- (n/(sum(w)^2)) * ((n*sum(w^2 * x^2) - sum(w*x)^2) / (n*(n-1)))
  ## Cochrane (1977)
  } else if (se_method=="cochrane") {
    mean_W <- mean(w)
    se <- (n/((n-1)*sum(w)^2))*(sum((w*x - mean_W*mean_Xw)^2)
                                - 2*mean_Xw*sum((w-mean_W)*(w*x-mean_W*mean_Xw))
                                + mean_Xw^2*sum((w-mean_W)^2))
  ## As implemented in Hmisc
  } else if (se_method=="Hmisc") {
    se <- (sum(w * (x - mean_Xw)^2) / (sum(w) - 1)) / n
  }
  return(sqrt(se))
}

where x is your vector of interest, w is a vector of weights with equal length, se_method specifies which method to use and na.rm specifies whether to remove missing values before performing calculations.

I understand that this doesn't fully answer your questions, but it might still be helpful to you.

You have multiple recursive (in the mathematical relation sense, not the computer science sense) expressions, notably on momentum, mean and variance. There is some minor vectorisation that can be done. I have tested this suggested code for correctness according to your provided examples. I've shown it to pull apart and illustrate what can be vectorised and what can't. I've also found accumulate to be a pain, because it seems to ignore the identity and dtype parameters.

Lo, there is hope: despite this looking as ugly as heck, in my testing it offers a ~180% speedup.

Suggested code

Docstrings omitted for brevity.

from timeit import timeit

import pandas as pd
import seaborn as sns
import numpy as np
from matplotlib import pyplot as plt
from numpy.random import default_rng


def ewm_old(values, time_steps, weights, decay_rate=0.01, variance_epsilon=1e-6):
    momentum = 0.0
    mean = 0.0
    variance = 0.0

    means = np.empty_like(values)
    stds = np.empty_like(values)
    for i in range(values.shape[0]):
        retention = (1.0 - decay_rate) ** time_steps[i]

        new_momentum = momentum * retention + weights[i]

        means[i] = mean = \
            mean * (momentum / new_momentum) * retention + \
            values[i] * weights[i] / new_momentum

        deviation = values[i] - mean

        variance = \
            variance * (momentum / new_momentum) * retention + \
            np.square(deviation) * weights[i] / new_momentum

        momentum = new_momentum

        stds[i] = np.sqrt(variance + variance_epsilon)

    return means, stds


def ewm_new(
    values: np.ndarray,
    time_steps: np.ndarray,
    weights: np.ndarray,
    decay_rate: float = 0.01,
    variance_epsilon: float = 1e-6,
) -> tuple[
    np.ndarray,
    np.ndarray,
]:
    def make_momentum(momentum: float, i: int) -> float:
        return momentum * retention[i-1] + weights[i-1]
    def make_means(mean: float, i: int) -> float:
        return mean*coefficients[i-1] + values[i-1]*offsets[i-1]
    def make_variances(variance: float, i: int) -> float:
        return variance*coefficients[i-1] + deviations[i-1]*offsets[i-1]

    retention = (1 - decay_rate)**time_steps
    momenta = np.frompyfunc(
        make_momentum, nin=2, nout=1, identity=0,
    ).accumulate(np.arange(len(values))).astype(np.float64)

    factors = momenta * retention / weights
    coefficients = factors / (factors + 1)
    offsets = 1 / (1 + factors)
    means = np.frompyfunc(
        make_means, nin=2, nout=1, identity=0,
    ).accumulate(np.arange(1+len(values)))[1:].astype(np.float64)

    deviations = (values - means)**2
    variances = np.frompyfunc(
        make_variances, nin=2, nout=1, identity=0,
    ).accumulate(np.arange(1+len(values)))[1:].astype(np.float64)

    stds = np.sqrt(variances + variance_epsilon)
    return means, stds


def test() -> None:
    rand = default_rng(seed=0)
    values, time_steps, weights = rand.random((3, 10))
    means, stds = ewm_new(values, time_steps, weights)

    assert np.allclose(
        means,
        np.array((
            0.63696169, 0.33792459, 0.09563421, 0.06077860, 0.28410923,
            0.38241663, 0.44761338, 0.51056811, 0.51504314, 0.56303342,
        )),
        rtol=0, atol=1e-6,
    )

    assert np.allclose(
        stds,
        np.array((
            0.00100000, 0.06149956, 0.05598960, 0.05115476, 0.29145246,
            0.34006516, 0.29895069, 0.28304818, 0.26340344, 0.27797243,
        )),
        rtol=0, atol=1e-6,
    )

    values     = np.array((1, 2, 1, 2, 1, 2, 1, 2, 1, 2), dtype=np.float64)
    time_steps = np.array((1, 1, 2, 2, 1, 1, 2, 2, 1, 1), dtype=np.float64)
    weights    = np.array((1, 1, 1, 1, 5, 5, 1, 1, 1, 1), dtype=np.float64)
    means, stds = ewm_new(values, time_steps, weights)

    assert np.allclose(
        means,
        (
            1.0000000,  1.50251256, 1.33219236, 1.50379090, 1.21925230,
            1.5028668,  1.46816320, 1.50314792, 1.47179905, 1.50307306,
        ),
        rtol=0, atol=1e-6,
    )

    assert np.allclose(
        stds,
        (
            0.00100000, 0.35266091, 0.34585958, 0.39006561, 0.30556194,
            0.38630063, 0.39249905, 0.40020053, 0.40503148, 0.41104515,
        ),
        rtol=0, atol=1e-6,
    )


def profile() -> None:
    rand = default_rng(seed=0)

    scale = np.round(10**np.linspace(0, 4, 150)).astype(int)
    times = []

    for n in scale:
        values, time_steps, weights = rand.random((3, n))

        for method in (ewm_old, ewm_new):
            t = timeit(lambda: method(values, time_steps, weights), number=1)
            times.append((method.__name__, n, t))

    df = pd.DataFrame(times, columns=('method', 'n', 'time'))
    sns.lineplot(data=df, x='n', y='time', hue='method')
    plt.show()


if __name__ == '__main__':
    test()
    profile()