The Fuzzy RD design can conceptualized as a local IV model (that is, an instrumental variables regression with weights that decline as observations move away from the cutoff). You need to instrument for the treated indicator with a dummy for being above the cutoff, while controlling for the running variable $Z$ and the interaction of above-the-cutoff dummy and $Z$. This can be found on page 958 of the 2nd edition of "Adult" Wooldridge. You don't have weights and you are missing these interactions in your two models.
Here's a simulation in Stata that demonstrates this equivalence. We start by installing two RD commands and making some fake data:
. clear
. /* install two commands that do fuzzy RD */
. capture net install rdrobust, from(http://www-personal.umich.edu/~cattaneo/rdrobust/stata) replace
. capture ssc install rd
. /* Generate Fake Data and Weights */
. mat c=(1,.5\.5,1)
. set seed 10011979
. drawnorm e pretest, n(1000) corr(c) clear
(obs 1000)
. gen z=pretest-0 // z is running variable
. gen above=z>0 // above is above-the-cutoff indicator
. gen treated=cond(uniform()<.8,above,1-above) // treated indicator
. gen y=z-z^3+treated+e // define outcome y
. gen w=max(0,1-abs(z)) // define triangle kernel weight
Here's the IV estimate. Note how you can do the interaction on the fly by using the factor variable notation. I am not using powers of $Z$ in my model, just a simple linear term:
. /* IV Version */
. ivregress 2sls y (treated=i.above) z c.z#i.above [pw=w]
(sum of wgt is 3.8525e+02)
Instrumental variables (2SLS) regression Number of obs = 703
Wald chi2(3) = 464.46
Prob > chi2 = 0.0000
R-squared = 0.3792
Root MSE = .9368
------------------------------------------------------------------------------
| Robust
y | Coef. Std. Err. z P>|z| [95% Conf. Interval]
-------------+----------------------------------------------------------------
treated | 1.641632 .3639524 4.51 0.000 .9282981 2.354965
z | .9815473 .2558127 3.84 0.000 .4801637 1.482931
|
above#c.z |
1 | -.4895004 .3488331 -1.40 0.161 -1.173201 .1942
|
_cons | -.2436488 .1800297 -1.35 0.176 -.5965006 .109203
------------------------------------------------------------------------------
Instrumented: treated
Instruments: z 1.above#c.z 1.above
There are two user-written commands that estimate fuzzy RD models:
. /* FRD Versions */
. rd y treated z, bw(1) z0(0) kernel(triangle)
Three variables specified; jump in treatment
at Z=0 will be estimated. Local Wald Estimate
is the ratio of jump in outcome to jump in treatment.
Assignment variable Z is z
Treatment variable X_T is treated
Outcome variable y is y
Estimating for bandwidth 1
Estimating for bandwidth .5
Estimating for bandwidth 2
------------------------------------------------------------------------------
y | Coef. Std. Err. z P>|z| [95% Conf. Interval]
-------------+----------------------------------------------------------------
numer | .7219002 .1669301 4.32 0.000 .3947233 1.049077
denom | .4397455 .0724289 6.07 0.000 .2977876 .5817035
lwald | 1.641632 .3642115 4.51 0.000 .9277902 2.355473
numer50 | .6432259 .240981 2.67 0.008 .1709118 1.11554
denom50 | .3158332 .1027168 3.07 0.002 .1145121 .5171544
lwald50 | 2.0366 .7731264 2.63 0.008 .5212996 3.5519
numer200 | 1.518334 .1283546 11.83 0.000 1.266764 1.769905
denom200 | .4920938 .0529886 9.29 0.000 .3882381 .5959496
lwald200 | 3.085456 .3284832 9.39 0.000 2.441641 3.729272
------------------------------------------------------------------------------
The first lwald coefficient is the FRD treatment effect. Here's another command that does FRD:
. rdrobust y z, fuzzy(treated) kernel(triangular) h(1) bwselect(IK)
Preparing data.
Computing variance-covariance matrix.
Computing RD estimates.
Estimation completed.
Sharp RD estimates using local polynomial regression.
Cutoff c = 0 | Left of c Right of c Number of obs = 1000
----------------------+---------------------- NN matches = 3
Number of obs | 373 330 BW type = Manual
Order loc. poly. (p) | 1 1 Kernel type = Triangular
Order bias (q) | 2 2
BW loc. poly. (h) | 1.000 1.000
BW bias (b) | 1.000 1.000
rho (h/b) | 1.000 1.000
Structural Estimates. Outcome: y. Running variable: z. Instrument: treated.
--------------------------------------------------------------------------------------
Method | Coef. Std. Err. z P>|z| [95% Conf. Interval]
----------------------+---------------------------------------------------------------
Conventional | 1.6416 .36912 4.4474 0.000 .918162 2.3651
Robust | - - 2.9955 0.003 .584165 2.79549
--------------------------------------------------------------------------------------
First-Stage Estimates. Outcome: treated. Running variable: z.
--------------------------------------------------------------------------------------
Method | Coef. Std. Err. z P>|z| [95% Conf. Interval]
----------------------+---------------------------------------------------------------
Conventional | .43975 .07148 6.1523 0.000 .299654 .579837
Robust | - - 2.7248 0.006 .083139 .509228
--------------------------------------------------------------------------------------
The conventional coefficient above is the FRD treatment effect. Both FRD estimates and their standard errors match the LWIV.
Now for you second question. Here I may be on shakier grounds since I am less familiar with the literature. I am assuming that you want to estimate a single model for men and women to get a single estimate of the effect. There are two options to accomplish this. One is to estimate the two models and re-weight the estimates. It seems prudent not to weight by the overall gender-specific sample size, nor by the treated sample size. Personally, I like to make the weights proportional to the number of units within some range of the discontinuity for each group to make sure that the observations too far from the cutoff don't matter in determining the weights. You can use the bandwidth for that. Because the estimates from the two individual discontinuities will be independent, once you have variances for each estimate it is easy to get a variance for the combined estimate since the covariance is zero.
The other option is to re-center all the observations by gender, pool them, and then apply an estimator for a single discontinuity with the running variable now a relative one rather than an absolute one. The resulting estimate implicitly weights the various discontinuity estimates by the number of observations at the discontinuity in each case.
I think I prefer the former approach, since it allows the bandwidth to vary by gender and for any heterogeneity in the treatment effect to emerge.