I'm not sure I understand the model itself, but as far as I can see, this is not an issue with the model.

The problem is that the function you pass to boot is not returning anything. There are two reasons for this: First, there is no element in the summary of a glm that is called estimate, so summary(stage_2)$estimate is NULL. Second, even if the first problem didn't exist, there is no coefficient that has the name "x2". This is because x2 is an ordered factor, which is internally turned into dummies.

You can adjust your code as follows:

glm_2siv_coef <- function(data,indices){
  d <- data[indices,]
  stage_1 <- polr(x2 ~ x + x3, Hess=TRUE, data=d)
  stage_2 <- glm(y ~ x + x2 + x3 + resids(stage_1),
                 family = "quasibinomial", data=d)
  return(stage_2$coefficients)
}
boot.results <- boot(data=df1, statistic=glm_2siv_coef, R=100)
boot.results

You will need to decide for yourself which estimates are of interest to you, but this gives you bootstrap SEs for all your coefficients. Also note that I increased R to 100, because a single bootstrap sample doesn't make any sense. I kept it low for illustrative purposes, so you should consider increasing it.

Your variable x2 is an ordinal variable, so when you call it in the glm, you will have n-1 coefficients, because the default contrast for ordinal dependent variable is contrast.poly, and it transform it into its quadratic terms (see more here):

glm(y ~ x + x2 + x3 ,data=df1)

Call:  glm(formula = y ~ x + x2 + x3, data = df1)

Coefficients:
(Intercept)            x         x2.L         x2.Q         x2.C         x2^4  
 -5.561e-16    1.000e-01    1.156e-17   -1.352e-18   -2.349e-17   -2.255e-17  
         x3  
 -5.294e-18 

The code you got from the link is a summary on another R object, but for glm, if you only need the coefficients, there's no harm in returning all the coefficients:

glm_2siv_coef <- function(data,indices){
  d <- data[indices,]
  stage_1 <- polr(x2 ~ x + x3, Hess=TRUE, data=d)
  stage_2 <- glm(y ~ x + x2 + x3 + resids(stage_1),
                 family = "quasibinomial", data=d)
  coefficients(stage_2)
}
boot.results <- boot(data=df1, statistic=glm_2siv_coef, R=10)

You can look the the bootstrap results like this:

mat = boot.results$t
colnames(mat) = names(glm_2siv_coef(df1,1:nrow(df1)))

head(mat)
     (Intercept)         x         x2.L        x2.Q         x2.C          x2^4
[1,]   -2.298085 0.4555501 -0.019243084 -0.01222671  0.003672487  0.0009351758
[2,]   -2.287667 0.4527313  0.024065483 -0.01124712  0.005012817  0.0028092668
[3,]   -2.297377 0.4543327  0.042194399 -0.02041060  0.007764768 -0.0025922548
[4,]   -2.301298 0.4565478  0.009459445 -0.01958682  0.000313143 -0.0058883027
[5,]   -2.296350 0.4544981 -0.010628593 -0.01141557 -0.005545030 -0.0030184318
[6,]   -2.297909 0.4537336  0.004486517 -0.01576363  0.016248869 -0.0005906523
               x3 resids(stage_1)
[1,] 9.114145e-04    0.0015985900
[2,] 1.187623e-03   -0.0020883363
[3,] 1.442589e-03   -0.0027610646
[4,] 9.258333e-05   -0.0021435513
[5,] 1.963208e-03    0.0005474028
[6,] 2.218785e-03    0.0014221113

The x2 related ones will be:

mat[,2:5]

This answer is based on this link, answer by @jay.sf.

set.seed(2)
sandbox$Group <- as.factor(sandbox$Group)
reduced.form <- polr(Group ~ z + random_variable + year, data=sandbox)
consistent.glm <- glm(y ~ Group + resids(reduced.form) + random_variable + year, family="quasibinomial", data=sandbox)

FUN <- function(x) {
  reduced.form <-  polr(Group ~ z + random_variable + year, data=x)
  fit <- glm(y ~ Group + resids(reduced.form) + random_variable + year,family="quasibinomial", data=x)
  fit$coefficients
}

set.seed(42)
R <- 200
bs <- t(replicate(R, FUN(sandbox[sample(nrow(sandbox), nrow(sandbox), replace=T), ])))

# To scrape out a summary, the matrixStats package is most convenient.

library(matrixStats)
b <- consistent.glm$coefficients
SE <- colSds(bs)
z <- b/SE
p <- 2 * pt(-abs(z), df = Inf)
ci <- colQuantiles(bs, probs=c(.025, .975))
res <- signif(cbind(b, SE, z, p, ci), 4)
res

DATA

a    <- 2    # structural parameter of interest
b    <- 1    # strength of instrument
rho  <- 0.5  # degree of endogeneity

N    <- 1000
z    <- rnorm(N)
res1 <- rnorm(N)
res2 <- res1*rho + sqrt(1-rho*rho)*rnorm(N)
x    <- z*b + res1
ys   <- x*a + res2
d    <- (ys>0) #dummy variable
y    <- (10-(d*ys))/10
random_variable <- rnorm(100, mean = 0, sd = 1)

library(data.table)
DT_1 <- data.frame(y,x,z, random_variable)
DT_2 <- structure(list(ID = c(1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 
                              13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 
                              29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 
                              45, 46, 47, 48, 49, 50), year = c(1995, 1995, 1995, 1995, 1995, 
                                                                1995, 1995, 1995, 1995, 1995, 2000, 2000, 2000, 2000, 2000, 2000, 
                                                                2000, 2000, 2000, 2000, 2005, 2005, 2005, 2005, 2005, 2005, 2005, 
                                                                2005, 2005, 2005, 2010, 2010, 2010, 2010, 2010, 2010, 2010, 2010, 
                                                                2010, 2010, 2015, 2015, 2015, 2015, 2015, 2015, 2015, 2015, 2015, 
                                                                2015), Group = c("A", "A", "A", "A", "B", "B", "B", "B", "C", 
                                                                                 "C", "A", "A", "A", "A", "B", "B", "B", "B", "C", "C", "A", "A", 
                                                                                 "A", "A", "B", "B", "B", "B", "C", "C", "A", "A", "A", "A", "B", 
                                                                                 "B", "B", "B", "C", "C", "A", "A", "A", "A", "B", "B", "B", "B", 
                                                                                 "C", "C"), event = c(1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 
                                                                                                      1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 
                                                                                                      1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1), win_or_lose = c(-1, 
                                                                                                                                                                    -1, -1, -1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
                                                                                                                                                                    0, 0, 0, 0, 0, 0, 0, 0, 0, -1, -1, -1, -1, 1, 1, 1, 1, 0, 0, 
                                                                                                                                                                    -1, -1, -1, -1, 1, 1, 1, 1, 0, 0)), row.names = c(NA, -50L), class = c("tbl_df", 
                                                                                                                                                                                                                                           "tbl", "data.frame"))
DT_1 <- setDT(DT_1)
DT_2 <- setDT(DT_2)
DT_2 <- rbind(DT_2 , DT_2 [rep(1:50, 19), ])
sandbox <- cbind(DT_1, DT_2)
sandbox <- setDT(sandbox)[y<0, y:=0]