How to write a block-bootstrap command in Stata?

Question

I am trying to write a code for block bootstrap in Stata for a panel data with time-series structure, but I somehow cannot save the results from it, it does not produce any results even when I am trying to store them into a matrix. If I ran code once in Stata - it works well, selects time-series blocks for each ID and calculated two-stage regression. Can someone please help me with that?

program define block_bootstrap, rclass
          // Save the original dataset state
        preserve             
        // Initialize matrix to store results
         if "`c(matrix_list)'" == "" | strpos("results", "`c(matrix_list)'") == 0 {matrix results = J(1, 2, .)  }
        // Generate blocks for all IDs at once
        sort ID year_month
        bysort ID (year_month): gen time_order = _n
        gen block = ceil(time_order / 5)  // Block size of 5 (adjust as needed)
    
        // Resample blocks with replacement
        foreach b in `block' {
            preserve
            tempfile boot_data
            keep if block == `b'
            xtoprobit endog_var instrument1, vce(robust)
    
            matrix coef = e(b)
            scalar cut1 = coef[1, 2]
            scalar cut2 = coef[1, 3]
            di "Cut1: " cut1
            di "Cut2: " cut2
            // Generate predictions and residuals
            predict xb_emp, xb
            gen predicted_category = .
            replace predicted_category = 1 if xb_emp <= cut1
            replace predicted_category = 2 if xb_emp > cut1 & xb_emp <= cut2
            replace predicted_category = 3 if xb_emp > cut2
            gen resid= endog_var - predicted_category
    
            // Second stage: Ordered Probit
            xtoprobit dependent_var endog_var resid, vce(robust)
            matrix coef2 = e(b)
            matrix V = e(V)
            scalar beta_OT5 = coef2[1, 1]  // 
            scalar se_OT5 = sqrt(V[1, 1])  // 
            return scalar beta_OT5
            return scalar se_OT5
            matrix results = results \ (beta_OT5, se_OT5)
            restore
                }
        return scalar beta = (results[., 1])  
        return scalar se = (results[., 2])    
        restore
    end
    
    bootstrap beta =beta, se =se , reps(100) cluster(ID): block_bootstrap