I am trying to calculate bootstrap standard errors for many regressions. I am applying Two-Sample Instrumental Variable method, so I use two samples and had to write my own bootstrap program. I've read through https://stats.idre.ucla.edu/stata/fa...strap-program/ and constructed my code accordingly. I am able to run my bootstrap program for a single regression. However, I need to do it for various regressions and want to replicate more than one regressions at the same time. Following code is the working one:
Code:
local varlist ""avearn" "incom" "labinc" "hour_wage" "haus" "eqinc""
***Here I create matrices to store my regression results
forval v=1/6{
        local w : word `v' of `varlist'

        forval z=1/2{
        local a=1+(`z'-1)*3
        qui reg cimp_`w' fimp_`w' if cfk090==`z', robust
            matrix f_`w'_`z'=_b[fimp_`w']
            
        qui reg cimp_`w' fimp_`w'  if cfk090==`z' & cfk070<=35, robust 
            matrix f_`w'_`z'35=_b[fimp_`w']
            
        qui reg cimp_`w' mimp_`w'  if cfk090==`z', robust
            matrix m_`w'_`z'=_b[mimp_`w']
            
        qui reg cimp_`w' mimp_`w'  if cfk090==`z' & cfk070<=35, robust 
            matrix f_`w'_`z'35=_b[mimp_`w']
        }
    }

****Bootstrap program
    program define myboot, rclass
    use cleanhaus.dta, replace
    bsample
    reg cavearn i.cfe030 b8.cageor i.cfi130 b2011.fb010 if cfk090==1 & celig==1, robust
    matrix male_coef_ave=e(b)
    use modul2011haus.dta, replace
    bsample
    gen fimp_avearn=. if cmf050!=. 
    forval p=0/6{
        forval r=1/9{
        di "`p'__`r'"
        local pp=`p'+1
        local rr=`r'+16

        replace fimp_avearn= male_coef_ave[1,39]+male_coef_ave[1,`pp']+male_coef_ave[1,`rr'] if cmf070==`p' & cmf110==`r'
        }
    }
    simulate f_avearn_1=_b[fimp_avearn], reps(100) seed(12345): myboot
    return scalar coef=_b[fimp_avearn]
    end

******Simulations
simulate f_avearn_1, reps(100) seed(12345): myboot
bstat, stat(f_avearn_1) n(13000)
Below is the one that I've tried to acquire bootstrap standard errors for coefficients from multiple regressions.
Code:
program define myboot, rclass

 use cleanhaus.dta, replace
 bsample
 
    reg clabinc i.cfe030 b8.cageor i.cfi130 b2011.fb010 if cfk090==1 & celig==1, robust
    matrix male_coef=e(b)
    reg clabinc i.cfe030 b8.cageor i.cfi130 b2011.fb010 if cfk090==2 & celig==1, robust
    matrix fem_coef=e(b)
    reg cavearn i.cfe030 b8.cageor i.cfi130 b2011.fb010 if cfk090==1 & celig==1, robust
    matrix male_coef_ave=e(b)
    reg cavearn i.cfe030 b8.cageor i.cfi130 b2011.fb010 if cfk090==2 & celig==1, robust
    matrix fem_coef_ave=e(b)
    reg cincom i.cfe030 b8.cageor i.cfi130 b2011.fb010 if cfk090==1 & celig==1, robust
    matrix male_coef_inc=e(b)
    reg cincom i.cfe030 b8.cageor i.cfi130 b2011.fb010 if cfk090==2 & celig==1, robust
    matrix fem_coef_inc=e(b)
    reg chour_wage i.cfe030 b8.cageor i.cfi130 b2011.fb010 if cfk090==1 & celig==1, robust
    matrix male_coef_hour=e(b)
    reg chour_wage i.cfe030 b8.cageor i.cfi130 b2011.fb010 if cfk090==2 & celig==1, robust
    matrix fem_coef_hour=e(b)
    reg chg110 i.cfe030 b8.cageor i.cfi130 b2011.fb010 if cfk090==1, robust
    matrix male_coef_haus=e(b)
    reg chg110 i.cfe030 b8.cageor i.cfi130 b2011.fb010 if cfk090==2, robust
    matrix fem_coef_haus=e(b)
    reg ceqinc i.cfe030 b8.cageor i.cfi130 b2011.fb010 if cfk090==1, robust
    matrix male_coef_eqinc=e(b)
    reg ceqinc i.cfe030 b8.cageor i.cfi130 b2011.fb010 if cfk090==2, robust
    matrix fem_coef_eqinc=e(b)
    
 use modul2011haus.dta, replace
 bsample
    capture drop fimp_labinc mimp_labinc fimp_avearn mimp_avearn fimp_incom mimp_incom fimp_hour_wage mimp_hour_wage fimp_haus mimp_haus fimp_eqinc mimp_eqinc
    gen fimp_labinc=. if cmf050!=. 
    gen mimp_labinc=. if cmf060!=. 
    gen fimp_avearn=. if cmf050!=. 
    gen mimp_avearn=. if cmf060!=. 
    gen fimp_incom=. if cmf050!=. 
    gen mimp_incom=. if cmf060!=. 
    gen fimp_hour_wage=. if cmf050!=. 
    gen mimp_hour_wage=. if cmf060!=. 
    gen fimp_haus=. if cmf050!=. 
    gen mimp_haus=. if cmf060!=. 
    gen fimp_eqinc=. if cmf050!=. 
    gen mimp_eqinc=. if cmf060!=. 
    
    forval p=0/6{
        forval r=1/9{
        di "`p'__`r'"
        local pp=`p'+1
        local rr=`r'+16

        replace fimp_labinc= male_coef[1,39]+male_coef[1,`pp']+male_coef[1,`rr'] if cmf070==`p' & cmf110==`r'
        replace mimp_labinc= fem_coef[1,39]+fem_coef[1,`pp']+fem_coef[1,`rr'] if cmf080==`p' & cmf140==`r'

        replace fimp_avearn= male_coef_ave[1,39]+male_coef_ave[1,`pp']+male_coef_ave[1,`rr'] if cmf070==`p' & cmf110==`r'
        replace mimp_avearn= fem_coef_ave[1,39]+fem_coef_ave[1,`pp']+fem_coef_ave[1,`rr'] if cmf080==`p' & cmf140==`r'

        replace fimp_incom= male_coef_inc[1,39]+male_coef_inc[1,`pp']+male_coef_inc[1,`rr'] if cmf070==`p' & cmf110==`r'
        replace mimp_incom= fem_coef_inc[1,39]+fem_coef_inc[1,`pp']+fem_coef_inc[1,`rr'] if cmf080==`p' & cmf140==`r'

        replace fimp_hour_wage= male_coef_hour[1,39]+male_coef_hour[1,`pp']+male_coef_hour[1,`rr'] if cmf070==`p' & cmf110==`r'
        replace mimp_hour_wage= fem_coef_hour[1,39]+fem_coef_hour[1,`pp']+fem_coef_hour[1,`rr'] if cmf080==`p' & cmf140==`r'

        replace fimp_haus= male_coef_haus[1,39]+male_coef_haus[1,`pp']+male_coef_haus[1,`rr'] if cmf070==`p' & cmf110==`r'
        replace mimp_haus= fem_coef_haus[1,39]+fem_coef_haus[1,`pp']+fem_coef_haus[1,`rr'] if cmf080==`p' & cmf140==`r'
        
        replace fimp_eqinc= male_coef_eqinc[1,39]+male_coef_eqinc[1,`pp']+male_coef_eqinc[1,`rr'] if cmf070==`p' & cmf110==`r'
        replace mimp_eqinc= fem_coef_eqinc[1,39]+fem_coef_eqinc[1,`pp']+fem_coef_eqinc[1,`rr'] if cmf080==`p' & cmf140==`r'

        }
    }

    local varlist ""avearn" "incom" "labinc" "hour_wage" "haus" "eqinc""
    forval v=1/6{
        local w : word `v' of `varlist'

        forval z=1/2{
        local a=1+(`z'-1)*3
        qui reg cimp_`w' fimp_`w' if cfk090==`z', robust
            return scalar f_`w'_`z'=_b[fimp_`w']
            
        qui reg cimp_`w' fimp_`w'  if cfk090==`z' & cfk070<=35, robust 
            return scalar f_`w'_`z'35=_b[fimp_`w']
            
        qui reg cimp_`w' mimp_`w'  if cfk090==`z', robust
            return scalar m_`w'_`z'=_b[mimp_`w']
            
        qui reg cimp_`w' mimp_`w'  if cfk090==`z' & cfk070<=35, robust 
            return scalar f_`w'_`z'35=_b[mimp_`w']
        }
    }    
end
simulate f_avearn_1=_b[fimp_avearn], reps(100) seed(12345): myboot
I am getting the error
"[fimp_avearn] not found
error in expression: _b[fimp_avearn]"
I would also appreciate if you point out any redundancies.

Thank you very much in advance!