Generalised Residuals - Estimated value of the density at that point for observation $j$ $i$?

Question

I am trying to write R code for the generalised residuals as described in Vella (1993) and also mentioned by Wooldridge (2014). Wooldridge suggests to use these residuals in order to estimate a discrete (ordinal) endogenous explanatory variable with the Control Function (CF) / Two Stage Residual Inclusion (2SRI) approach (see Wooldridge, 2014 and Terza et al. 2008).

They seem pretty important, as it appears to be the only way to somewhat consistently estimate an ordinal endogenous explanatory variable in a non linear structural equation. Nevertheless they are not available in either R or Stata (EDIT: they are available in stata for oprobit, by using predict gr, score). As a result I'm trying to have a go myself.

The residual is described as follows:

Generalised Residuals

UPDATE: According to this post, $\pi$, is the numerator of the formula (3) below, which comes from Chiburis and Lokshin (2007), where $\Pi$ is the denominator.

$W_i$ are the independent variables

$z_i$ is the categorical variable

$\alpha$ is an unknown vector of parameters

Question

I started with an attempt to manually calculate the residuals, but I do not really understand what to do with $\pi_ji)$, which is estimated value of the density at that point. Could anyone tell me what that should be?

I think that is the only thing missing or unclear in the code below.. It would be great if someone could help me calculate $\pi$ from the polr output.

R code for the generalised residual with example data

library(sure) # for residual function and sample data sets
library(MASS) # for polr function
library(VGAM)
rm(df1, df2, df3)
df1 <- df1
df1$x2 <- df2$y
df1$x3 <- df2$x
df1$y <- df3$x/10
df1$x <- df1$x/0.8*df1$x

polr  <- polr(x2 ~ x + x3, method="probit", Hess=TRUE, data=df1) 

predicted_probabilities <- predict(polr, type="probs")

category_dummies <- model.matrix( ~ x2- 1, data=df1)

cum_probability <- polr$fitted.values
cum_probability <- t(apply(cum_probability, 1, cumsum)) # https://stackoverflow.com/questions/67301413/automatically-creating-a-cumulative-table-by-column?noredirect=1#comment118959242_67301413

residuals <- model.matrix( ~ x2 - 1, data=df1) # Just to make it match the df
for(i in seq_len(nrow(df1))){
 for(j in seq_len(ncol(df1))){
  residuals[i,j] <- category_dummies[i,j] * cum_probability[i,j] * ( 1 / predicted_probabilities[i,j]) * ( 1 / ( 1 - predicted_probabilities[i,j])) * ( predicted_probabilities[i,j] - predicted_probabilities[i,j] )
 }}

In Stata

oprobit x2 c.x c.x3

Iteration 0:   log likelihood = -2052.0298  
Iteration 1:   log likelihood = -2039.4155  
Iteration 2:   log likelihood = -2039.4084  
Iteration 3:   log likelihood = -2039.4084  

Ordered probit regression                       Number of obs     =      2,000
                                                LR chi2(2)        =      25.24
                                                Prob > chi2       =     0.0000
Log likelihood = -2039.4084                     Pseudo R2         =     0.0062

------------------------------------------------------------------------------
          x2 |      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
           x |   -.001509   .0014161    -1.07   0.287    -.0042845    .0012665
          x3 |   .0884805   .0180269     4.91   0.000     .0531485    .1238126
-------------+----------------------------------------------------------------
       /cut1 |  -1.588152   .1101114                     -1.803967   -1.372338
       /cut2 |  -.7587626   .0933795                     -.9417829   -.5757422
       /cut3 |   1.227218   .0948042                      1.041405    1.413031
       /cut4 |   2.057595   .1074906                      1.846918    2.268273
------------------------------------------------------------------------------

predict gr, score

References

Chiburis, R., & Lokshin, M. (2007). Maximum Likelihood and Two-Step Estimation of an Ordered-Probit Selection Model. The Stata Journal, 7(2), 167–182. https://doi.org/10.1177/1536867X0700700202

Terza J. V, Basu A., Rathouz, P.J. (2008) Two-stage residual inclusion estimation: addressing endogeneity in health econometric modeling. Journal of health economics 27 (3), 531-543, 2008

Vella, F. (1993). A Simple Estimator for Simultaneous Models with Censored Endogenous Regressors. International Economic Review, 34(2), 441-457. doi:10.2307/2526924

Wooldridge,JM. (2014) Quasi-maximum likelihood estimation and testing for nonlinear models with endogenous explanatory variables. Journal of Econometrics Volume 182, Issue 1, September 2014, Pages 226-234