Returning to answer my own question, one way to estimate the variance is to use M-estimation (or estimating equations as stated by Noah). This method relies on the use of parametric models for $E[Y|A,W,V]$ and $\Pr(A | W,V)$. Below is a description of how this can be done and an example.
M-estimation
For an intro to M-estimation, see Stefanski & Boos (2002) or Cole et al. (2022). Essentially, we are going to simultaneously solve a series of estimating equations. Then we use the sandwich variance to estimate the variance for $\beta$.
The stacked estimating equations are
$$ \sum_{i=1}^{n} \psi(O_i; \theta) = \sum_{i=1}^{n}
\begin{bmatrix}
\left(A_i - \text{expit}(X_i^T \alpha)\right)X_i \\
\left(Y_i - Z_i^T \gamma\right)Z_i \\
\left[(\tilde{Y}_i^1 - \tilde{Y}_i^0) - V_i^T \beta\right]Z_i \\
\end{bmatrix} = 0
$$
where $O_i = (Y_i, A_i, W_i, V_i)$, $\theta = (\alpha, \gamma, \beta)$, $X_i$ is the design matrix for the propensity score model, $Z_i$ is the design matrix for the outcome model, and
$$\tilde{Y}_i^a = \frac{Y_i I(A_i=a)}{A_i \text{expit}(X_i^T \alpha) + (1-A_i)(1-\text{expit}(X_i^T \alpha))} - \frac{{Z_i^a}^T (\text{expit}(X_i^T \alpha) - A_i)}{\text{expit}(X_i^T \alpha)}$$
where $Z_i^a$ indicating the design matrix with $A_i=a$. Note that $\tilde{Y}_i^a$ are the pseudo-potential-outcomes under $a$ mentioned in step 3 of the question.
The first estimating equation is the propensity score model (i.e., logistic regression), the second is the outcome model (i.e., linear regression), and the final is the CATE above. To solve the estimating equations, we use a root-finding procedure to find the values of $\hat{\theta}$ that result in the sum of the estimating functions being zero (or nearly so).
After, the empirical sandwich variance estimator is used to estimate the full covariance matrix (and the variance for $\beta$). The sandwich is
$$V(O_i, \hat{\theta}) = B(O_i, \hat{\theta})^{-1} M(O_i, \hat{\theta}) \left\{ B(O_i, \hat{\theta})^{-1} \right\}^T$$
where
$$B(O_i, \hat{\theta}) = n^{-1} \sum_i -\psi'(O_i, \hat{\theta})$$
with $\psi'$ indicating the matrix derivative and
$$M(O_i, \hat{\theta}) = n^{-1} \sum_i \psi(O_i, \hat{\theta}) \psi(O_i, \hat{\theta})^T$$
This variance estimator allows the uncertainty of $\hat{\beta}$ to depend on the uncertainty of $\hat{\alpha}$ and $\hat{\gamma}$.
Code
To show how this would be implemented programmatically, the following is a simple example using the Python M-estimation library, delicatessen (disclosure: I am the creator of delicatessen).
Initializing
import numpy as np
import pandas as pd
from delicatessen import MEstimator
from delicatessen.estimating_equations import ee_regression
from delicatessen.utilities import inverse_logit
Generating some generic data
np.random.seed(20230109)
n = 1000
d = pd.DataFrame()
d['W'] = np.random.normal(size=n)
d['V'] = np.random.binomial(n=1, p=0.5, size=n)
pr_a = logistic.cdf(-2 + 2*d['V'] + d['W'])
d['A'] = np.random.binomial(n=1, p=pr_a, size=n)
d['Y1'] = 3 + 1*1 + 1*1*d['V'] + d['V'] - 1.2*d['W'] + np.random.normal(size=n)
d['Y0'] = 3 + 0*1 + 0*1*d['V'] + d['V'] - 1.2*d['W'] + np.random.normal(size=n)
d['Y'] = np.where(d['A'] == 1, d['Y1'], d['Y0'])
d['intercept'] = 1
d['AV'] = d['A']*d['V']
Setting up design matrices as NumPy arrays for delicatessen
X = np.asarray(d[['intercept', 'V', 'W']])
a = np.asarray(d['A'])
Z = np.asarray(d[['intercept', 'A', 'AV', 'V', 'W']])
y = np.asarray(d['Y'])
V = np.asarray(d[['intercept', 'V']])
da = d.copy()
da['A'] = 1
da['AV'] = da['A']*d['V']
Z1 = np.asarray(da[['intercept', 'A', 'AV', 'V', 'W']])
da['A'] = 0
da['AV'] = da['A']*d['V']
Z0 = np.asarray(da[['intercept', 'A', 'AV', 'V', 'W']])
Defining stacked estimating equations
def psi(theta):
# Separating theta into the constituent components
alpha = theta[0:3] # Prop score params
gamma = theta[3:8] # Outcome params
beta = theta[8:] # CATE params
# Estimating equations for propensity score model
ee_ps = ee_regression(theta=alpha, # Model parameters
X=X, # ... X design matrix
y=a, # ... treatment
model='logistic') # ... with logistic model
# Estimating equations for outcome model
ee_out = ee_regression(theta=gamma, # Model parameters
X=Z, # ... Z design matrix
y=y, # ... observed outcome
model='linear') # ... with linear model
# Generating pseudo potential outcomes
pi_a = inverse_logit(np.dot(X, alpha)) # Propensity score
y_a1 = np.dot(Z1, gamma) # Predicted Y^1
y_a0 = np.dot(Z0, gamma) # Predicted Y^0
yt_a1 = y*a/pi_a - y_a1*(a-pi_a)/pi_a # Pseudo Y^1
yt_a0 = y*(1-a)/(1-pi_a) + y_a0*(a-pi_a)/(1-pi_a) # Pseudo Y^0
# Estimating CATE model
ee_cate = ee_regression(theta=beta, # Parameters
X=V, # ... with CATE design matrix
y=yt_a1 - yt_a0, # ... difference is pseudo outcomes
model='linear') # ... with linear model
# Returning stack of estimating equations
return np.vstack([ee_ps, ee_out, ee_cate])
Applying the M-estimator
estr = MEstimator(psi, init=[0, ]*10)
estr.estimate(solver='lm')
estr.theta[8:] # Point estimates
estr.variance[8:, 8:] # Covariance matrix
which results in $\hat{\beta}_0 = 1.00$ (variance: 0.020) and $\hat{\beta}_1 = 0.85$ (variance: 0.029).
More examples of M-estimators and an introduction to delicatessen is provided in Zivich et al. (2022).
Caveat
This approach may not work when machine learning is used to estimate the nuisance models. Instead, the influence function could be used to derive a variance estimator.
References
Cole, S. R., Edwards, J. K., Breskin, A., Rosin, S., Zivich, P. N., Shook-Sa, B. E., & Hudgens, M. G. (2022). Illustration of Two Fusion Designs and Estimators. American Journal of Epidemiology.
Stefanski, L. A., & Boos, D. D. (2002). The calculus of M-estimation. The American Statistician, 56(1), 29-38.
Zivich, P. N., Klose, M., Cole, S. R., Edwards, J. K., & Shook-Sa, B. E. (2022). Delicatessen: M-Estimation in Python. arXiv preprint arXiv:2203.11300.
