# semiparametric index model with heteroskedasticity

I'm trying to estimate a semiparametric binary response model with index heteroscedasticity in R. That is, I have a model defined with $$y_i = \mathbf{1}\{\beta_0 + \beta_1 x_{1i} + \beta_2 x_{2i} + \epsilon_{i} >0\}$$ where $$\mathbf{1}\{\cdot \}$$ is the indicator function, so $$y_i \in \{0,1\}$$, and $$\epsilon \sim G(0, \sigma(x_1)^2)$$, where $$G$$ is some distribution function with mean zero and variance $$\sigma(x_1)^2 = \exp(\delta x_1).$$ I'm interested in estimating $$\beta$$'s, and $$\delta$$ together with function $$G$$. So, you can think of the model as say a heteroskedastic logit model, where instead of fixing logit function I would like to additionally estimate $$G$$.

I was trying to use np package in R. The package contains implemented function for 'Klein & Spady' single index semiparmetric estimator (npindexbw, npindex). Here is an example:

require(np)
set.seed(12345)
n <- 1000 # no.
const=runif(n,1,1)
x1 <- runif(n,-1,1) # predictor 1
x2 <- runif(n,-1,1) # " 2
e1 <- rnorm(n,0,1) # normal error
e2 <- (0.5+ 0.5*(x1))*e1 # heteroskedastic error
y <- ifelse(0.5 + 0.5*x1 -0.5*x2 - e2 >0, 1, 0) #outcome

#Estimation:

However, npindexbw doesn't allow for the explicit specification of the heteroskedasticity, so cannot estimate separately $$\beta$$'s and $$\delta$$. Any suggestion on how to estimate the model of the type described above?