3 question for the smartest on using and interpreting Oster's -psacalc- from SSC

[Peter Meier]

Dear smartest listers,

I used to estimate

reghdfe y x1##x2, absorb(a##b)

yielding an estimated coefficient on x1#x2 of 0.15.

Since I am concerned that selection into the x2 treatment may be endogenous to some unobservable, I next added a control to estimate

reghdfe y x1##x2 x1##x3, absorb(a##b)

which, rather than reducing the coefficient of interest has slightly increased it to 0.20.

Now I understand from Oster (2019) that before I can satisfied with the coefficient being relatively stable and if anything larger when controlling for the other interaction,
I should take into account also how much the inclusion increases my R2 and the variance of the controls,
which I understand can be implemented with the user-written -psacalc- available from SSC.

As psacalc seems not to work after -reghdfe- but after -areg-, and -areg- seems not to like interactions, I have hence written the following lines:

gen x1_x2 = x1*x2
gen x1_x3 = x1*x3
egen abfe = group(a b)
areg y x1_x2 x1 x2 x1_x3 x1 x3, absorb(abfe)
psacalc beta x1_x2

Now I face the following issues:

(1) It tells me the controlled coefficient is 0.20 (correct) and the uncontrolled one is 0.70. When I check I find that 0.70 is what I would get by regressing y only on x1_x2 without the main terms, which I would never have considered doing. So in my view it benchmarks my controlled model against the wrong "uncontrolled" one, but I am not sure how to tell -psacalc- what to use as benchmark. I tried with

psacalc beta x1_x2, mcontrols(x1 x2 abfe)

to clarify that I want the main effects and the fixed effects always in regardless of the possible selection on unobservables,
but then I am just being told that "option mcontrols() not allowed"
Does anyone see where I am getting this wrong and how I can get psacalc to compare my controlled model with both x1_x2 and x1_x3 and all its main terms to the uncontrolled model which misses x1_x3 and its main terms but does still contain the main effects as well as the fixed effects?

(2) The command output then says Beta = -350, "Alt. sol. 1" = 0.8. Should I understand that it suggests the true beta to be -350 which would not only have the opposite sign of any coefficients I have obtained so far but also be about 50 times as large in absolute magnitude which seems totally weird to me. Or am I misreading this?

(3) Finally, next to its alternative solution of 0.8 it says "Bias changes direction: Yes". But why so when my controlled, my claimed and true uncontrolled, and the here proposed alternative solution are all positive? On that one I have already seen an earlier post by someone unsure how to read this "changes direction" but no response yet, so I hope someone can help this time?

Thanks so much!
PM

[George Ford]

areg doesn't work right, but xtreg does.

clear all

sysuse auto, clear
drop if mi(rep78)

eststo e1: reg weight c.foreign#c.mpg 
eststo e2: areg weight c.foreign#c.mpg , absorb(rep78)
eststo e3: areg weight c.foreign##c.mpg length, absorb(rep78)
esttab e1 e2 e3
psacalc beta c.foreign#c.mpg 

xtset rep78
eststo d1: reg weight c.foreign#c.mpg 
eststo d2: xtreg weight c.foreign#c.mpg , fe
eststo d3: xtreg weight c.foreign##c.mpg length, fe
esttab d1 d2 d3
psacalc beta c.foreign#c.mpg

[Peter Meier]

Dear George,

thank you so much! Indeed, as your example nicely demonstrates, psacalc seems to compare estimates 3 to estimates 2, as I want, when I precede it with xtreg, whereas it wrongly compares them to estimates 1 when I precede it with areg. This was not clear to me as the psacalc help file had mentioned both areg and xtreg as allowed estimation commands, very helpful that you discovered this!

While it has an option to explicitly specify the "controlled" model, it has no option to explicitly specify the "uncontrolled" one, so it is not transparent to me where it takes its benchmark from.
But I seem to get the desired comparison of d3 vs d2 also when not explicitly running d2, so it seems that rather than comparing the last 2 models run it compares the last run against the last run minus that stripped of all controls not specified in the "mcontrol" option (which works, while "mcontrols" as it is called in the help file produces an error).

As in all models I want to have not only the foreign*mpg interaction but also the main terms, I now run the following (clumsily including all interactions manually to make sure none is forgotten):

 
 clear all sysuse auto, clear drop if mi(rep78)  
 gen foreign_mpg = foreign*mpg  xtset rep78 eststo d1: reg weight foreign_mpg foreign mpg  *eststo d2: xtreg weight foreign_mpg foreign mpg , fe eststo d3: xtreg weight foreign_mpg foreign mpg length, fe esttab d1 /*d2*/ d3, s(N R2) psacalc beta c.foreign#c.mpg, mcontrol(foreign mpg)

That does indeed seem to compare the right coefficients and the right R2.

I still have 2 remaining concerns, not sure if you or anyone else here can help with them:

1. -xtreg- seems to produce way smaller R2 than -reghdfe-. As I understand the R2 to be a key input to -psacalc- I gather this might matter but am not sure what to do about it?
2. In the above and other examples, the "beta" reported by -psacalc- has a sign different from any of the coefficients I estimate with or without controls, while the "Alt. sol. 1" does not. And I do not fully understand from the psacalc help file which one I should report and rely on? Maybe only the providers of that ado can help with this, but maybe one of you understands it better than I do as well?

Thanks so much,
PM

[George Ford]

You want to use Beta, not Alt. Sol. 1. Oster's calculations can have multiple solutions, so Stata chooses the best one.

The code could use an update for multiple fixed effects and to remedy some odd behaviors. It does not look like areg is incorporating the fixed effects properly. It does not work properly if you have multiple fe in xtreg. You can get it to work with mcontrol() but have to use xi: and then include all the created fe.

I may take a shot at it at some point, but there are coders that could do a better job.

[George Ford]

Try this. It could be cleaned up a bit but it works.

Here's how it goes. You specify a constrained model (give it a name, here "e1"), and a unconstrained model ("e2"). You can name them what you want but I think the constrained model must be first.

eststo e1: reghdfe y c.post#c.treated , absorb(id t)
eststo e2: reghdfe y c.post#c.treated x1 , absorb(id t)
reghdfe_psacalc beta c.post#c.treated , models(e1 e2)

[George Ford]

reghdfe_psacalc.do

[George Ford]

clear all

set obs 20
g id = _n
g fe = rchi2(2)
expand 10
bys id: g t = _n

g treated = id <= 5
g post = t > 5

g x1 = rgamma(5,1)

g y = fe + 0.5*x1 + 2*post*treated + rnormal()

xtset id t
xtreg y c.post#c.treated i.t , fe
xi: xtreg y c.post#c.treated x1 i.t, fe
psacalc beta c.post#c.treated , mcontrol(_It_2 _It_3 _It_4 _It_5 _It_6 _It_7 _It_8 _It_9 _It_10)

eststo e1: reghdfe y c.post#c.treated , absorb(id t)
eststo e2: reghdfe y c.post#c.treated x1 , absorb(id t)
reghdfe_psacalc beta c.post#c.treated , models(e1 e2)

[George Ford]

Note: I used Claude.ai to rewrite psacalc. It took some back and forth, but it eventually worked. Whether it handles more complex issues (delta option, eg) is untested. It could be coded to extract the covariates, but giving it the two models seemed a better place to start and allows the avoidance of mcontrol(), though mcontrol() option is still there, though I'm not sure if that works correctly.

[Peter Meier]

Thank you! Major achievement if one can run it also after reghdfe and if one can explicitly specify which 2 models to compare to each other.

Unfortunately on my end Stata on one computer keeps saying "command reghdfe_psacalc is unrecognized", even though I saved it into the same personal folder where I had also placed the psacalc and not gotten this error and even though I tried saving it both as "do" and as "ado" files. On the other computer on which I tried, I instead got the error ""/ is not a valid command name" when trying to run it.

My local AI program pointed out that an ado file should start with "program define ..." rather than just with "program ...". Just to make sure this is not the issue I added "define" but that did not solve it. Googling it also suggested that Claude.ai may have added characters not recognizable by Stata and invisible to the human eye, but when I open the text in Notepad++ it says it is encoded as "UTF-8" which seems to be ok for Stata. Any idea on why I might be getting these errors when it seems to have run through on your end?

Thanks so much, PM

[Nick Cox]

My local AI program pointed out that an ado file should start with "program define ..." rather than just with "program ...".

This is quite wrong. The help for program alone shows that define is optional there. If memory serves me right, define was necessary long, long ago, but not for many years (versions).

As a detail on the side, the expression for the smartest is in my view better avoided. I take it that you don't want to put down or discourage anyone at all interested in reading or responding to your question.

I am not smart enough myself to follow exactly what you've done. Sorry. I'd recommend copying here the file in question.

[Peter Meier]

Thanks, Nick.

With my "smartest" framing I just wanted to express that I personally could not solve this and have a lot of respect for all those who can. If you think it might to the contrary discourage readers who like me cannot then I apologize.

Happy to post below the code which George thankfully posted above to ensure that the existing -psacalc- can use also inputs from -reghdfe- rather than just from -xtreg- and that one can explicitly specify which 2 models to compare:

capture program drop reghdfe_psacalc    


*! reghdfe_psacalc 1.0 27Sep2025
 
       
program reghdfe_psacalc, rclass

    version 11.2
    syntax anything [, models(string) rmax(real 1.0) delta(real 1.0) beta(real 0)]
    
    * Get treatment var and type of calculation
    
    tokenize `anything'
    local type "`1'"
    local treatment "`2'"
    
    if "`3'"!="" {
        di as err _n "Too many variables"
        exit 103    
    }
    
    * Parse models option
    if "`models'"=="" {
        di as err _n "models() option required: specify models(constrained_model full_model)"
        exit 198
    }
    
    tokenize `models'
    local constrained_model "`1'"
    local full_model "`2'"
    
    if "`constrained_model'"=="" | "`full_model'"=="" {
        di as err _n "models() must specify exactly two model names: models(constrained_model full_model)"
        exit 198
    }
    
    if "`3'"!="" {
        di as err _n "models() must specify exactly two model names"
        exit 198
    }
    
    tempvar error_hat treat_res touse treatment_var
    tempname beta_tilde r_tilde beta_o r_o sigma_yy sigma_xx t_x scbeta scdelta scrmax
    
    /*=========================================================================
                        1: Validate stored models and capture errors
    ===========================================================================*/
    
    * Check if models exist and are reghdfe
    cap estimates restore `constrained_model'
    if _rc {
        di as err _n "Constrained model `constrained_model' not found"
        exit 301
    }
    
    if "`e(cmd)'"!="reghdfe" {
        di as err _n "Constrained model must be estimated with reghdfe"
        exit 301
    }
    
    cap estimates restore `full_model'
    if _rc {
        di as err _n "Full model `full_model' not found"
        exit 301
    }
    
    if "`e(cmd)'"!="reghdfe" {
        di as err _n "Full model must be estimated with reghdfe"
        exit 301
    }
    
    * Start with full model (controlled)
    estimates restore `full_model'
    local names_full: colnames e(b)
    local depvar_full = e(depvar)
    scalar `beta_tilde' = _b[`treatment']
    scalar `r_tilde' = e(r2)
    gen `touse' = e(sample)
    
    * Get weights if regression had them
    if "`e(wtype)'"!="" {
        loc weight "`=e(wtype)' `=e(wexp)'"
    }
    
    * Switch to constrained model (uncontrolled)
    estimates restore `constrained_model'
    local names_constrained: colnames e(b)
    local depvar_constrained = e(depvar)
    
    * Validate models are compatible
    if "`depvar_full'"!="`depvar_constrained'" {
        di as err _n "Models must have the same dependent variable"
        exit 198
    }
    
    scalar `beta_o' = _b[`treatment']
    scalar `r_o' = e(r2)
    
    // Delta ignored if requesting delta and entered delta!=1
    if "`type'"=="delta" & `delta'!=1 di as error "User entered delta ignored"
    
    // Beta ignored if requesting beta and entered beta!=1
    if "`type'"=="beta" & `beta'!=0 di as error "User entered beta ignored"        
    
    // Does not have a constant in full model
    if strpos("`names_full'","_cons")==0  {
        di as error "Full model does not have a constant"
        exit 5001
    }
    
    // Does not have a constant in constrained model
    if strpos("`names_constrained'","_cons")==0  {
        di as error "Constrained model does not have a constant"
        exit 5001
    }
    
    // Treatment not in full model
    if strpos("`names_full'","`treatment'")==0 {
        di as err _n "Treatment variable `treatment' not found in full model"
        exit 5001
    }
    
    // Treatment not in constrained model
    if strpos("`names_constrained'","`treatment'")==0 {
        di as err _n "Treatment variable `treatment' not found in constrained model"
        exit 5001
    }
    
    // rmax out of bounds
    if `rmax'>1 {
        di as err _n "The maximum possible R-squared is 1"
        exit 5001
    }
    
    // rmax less than controlled R-squared
    if `rmax'<`r_tilde' {
        di as err _n "Maximum R-squared provided is less than controlled R-squared"
        exit 5001
    }
    
    // bad type
    if "`type'"~="beta" & "`type'"~="delta" {
        di as err _n "Invalid type: must be beta or delta"
        exit 5001
    }
    
    /*=========================================================================
                        2: Calculate variance components
    ===========================================================================*/
    
    * Get dependent variable variance
    qui sum `depvar_full' if `touse'
    scalar `sigma_yy' = r(Var)
    
    * Handle factor variables - create temporary treatment variable  
    * For factor variables like c.post#c.treated, we need to create the actual variable
    cap confirm variable `treatment'
    if _rc {
        * This is a factor variable expression - parse it
        * For c.post#c.treated, we want post*treated
        local factor_expr = "`treatment'"
        
        * Remove c. prefixes and replace # with *
        local factor_expr = subinstr("`factor_expr'", "c.", "", .)
        local factor_expr = subinstr("`factor_expr'", "#", "*", .)
        
        * Generate the interaction variable
        qui gen `treatment_var' = `factor_expr' if `touse'
    }
    else {
        * Treatment is a regular variable
        qui gen `treatment_var' = `treatment' if `touse'
    }

    * Get treatment variable variance after residualizing on constrained controls
    estimates restore `constrained_model'
    local constrained_vars: colnames e(b)
    local cons "_cons"
    local constrained_controls: list constrained_vars - cons
    local constrained_controls: list constrained_controls - treatment
    local constrained_absorb = e(absvars)
    
    * Residualize treatment on constrained controls
    if "`constrained_controls'"!="" {
        qui reghdfe `treatment_var' `constrained_controls' [`weight'] if `touse', absorb(`constrained_absorb') resid
    }
    else {
        qui reghdfe `treatment_var' [`weight'] if `touse', absorb(`constrained_absorb') resid
    }
    qui predict `treat_res' if `touse', resid
    qui sum `treat_res' if `touse'
    scalar `sigma_xx' = r(Var)
    
    * Get variance of treatment residualized on all controls (for t_x)
    estimates restore `full_model'
    local full_vars: colnames e(b)
    local full_controls: list full_vars - cons
    local full_controls: list full_controls - treatment
    local full_absorb = e(absvars)
    
    qui reghdfe `treatment_var' `full_controls' [`weight'] if `touse', absorb(`full_absorb') resid
    qui predict `error_hat' if `touse', resid
    qui sum `error_hat' if `touse'
    scalar `t_x' = r(Var)
    
    scalar `scbeta' = `beta'
    scalar `scdelta' = `delta'
    scalar `scrmax' = `rmax'
    
    /*=========================================================================
                        3: Solve for outputs in mata
    ===========================================================================*/

    mata {
        beta_o = st_numscalar("`beta_o'")    
        beta_tilde = st_numscalar("`beta_tilde'")
        r_o = st_numscalar("`r_o'")
        r_tilde = st_numscalar("`r_tilde'")
        sigma_yy = st_numscalar("`sigma_yy'")
        sigma_xx = st_numscalar("`sigma_xx'")
        t_x = st_numscalar("`t_x'")
        delta = st_numscalar("`scdelta'")
        beta = st_numscalar("`scbeta'")
        r_max = st_numscalar("`scrmax'")    
        
        bo_m_bt = beta_o - beta_tilde
        rt_m_ro_t_syy = (r_tilde - r_o)*sigma_yy
        rm_m_rt_t_syy = (r_max - r_tilde)*sigma_yy
        
        // Solve delta bound
        bt_m_b = beta_tilde - beta
        num1 = (bt_m_b)*rt_m_ro_t_syy*t_x
        num2 = (bt_m_b)*sigma_xx*t_x*(bo_m_bt)^2        
        num3 = 2*(bt_m_b)^2*(t_x*bo_m_bt*sigma_xx)
        num4 = ((bt_m_b)^3)*((t_x*sigma_xx-t_x^2))
        
        den1 = rm_m_rt_t_syy*bo_m_bt*sigma_xx
        den2 = bt_m_b*rm_m_rt_t_syy*(sigma_xx-t_x)
        den3 = (bt_m_b^2)*(t_x*bo_m_bt*sigma_xx)
        den4 = (bt_m_b^3)*(t_x*sigma_xx-t_x^2)
        
        num = num1+num2+num3+num4
        den = den1+den2+den3+den4
        
        ds = num/den
        st_numscalar("boundx", ds)
        
        // Solve for beta if delta = 1
        if (delta==1) {
            cap_theta = rm_m_rt_t_syy*(sigma_xx-t_x)-rt_m_ro_t_syy*t_x-sigma_xx*t_x*(bo_m_bt^2)
            d1_1 = 4*rm_m_rt_t_syy*(bo_m_bt^2)*(sigma_xx^2)*t_x
            d1_2 = -2*t_x*bo_m_bt*sigma_xx
            
            sol1 = (-1*cap_theta-sqrt((cap_theta)^2+d1_1))/(d1_2)
            sol2 = (-1*cap_theta+sqrt((cap_theta)^2+d1_1))/(d1_2)
            
            beta1 = beta_tilde - sol1
            beta2 = beta_tilde - sol2
            
            if ( (beta1-beta_tilde)^2 < (beta2-beta_tilde)^2) {
                betax = beta1
                altsol1 = beta2
            }
            else {
                betax = beta2
                altsol1 = beta1
            }
            
            // Change to altsol if bias direction differs
            if ( sign(betax-beta_tilde)!=sign(beta_tilde-beta_o) ) {
                solc = betax
                betax = altsol1
                altsol1 = solc
            }
            
            distx = (betax - beta_tilde)^2
            dist1 = (altsol1 - beta_tilde)^2
            
            st_numscalar("betax", betax)
            st_numscalar("altroot1", altsol1)
            st_numscalar("altroot2", .)
            st_numscalar("distx", distx)
            st_numscalar("dist1", dist1)
            st_numscalar("dist2", .)
        }
        else {
            // For delta != 1, just return the main solution for now
            // (simplified version - full cubic solution can be added later)
            st_numscalar("betax", beta_tilde)
            st_numscalar("altroot1", .)
            st_numscalar("altroot2", .)
            st_numscalar("distx", 0)
            st_numscalar("dist1", .)
            st_numscalar("dist2", .)
        }
    }
    
    loc boundx = boundx
        
    if `delta'==1 {
        loc j=2
        foreach x in betax altroot1 distx dist1 {
            loc `x' = `x'
        }        
    }
    
    else if `delta'~=1 {
        foreach x in betax altroot1 altroot2 distx dist1 dist2 {
            loc `x' = `x'
        }
        if `altroot1'==. & `altroot2'==. loc j=1
        else if (`altroot1'==. & `altroot2'!=. | `altroot2'!=. & `altroot2'==.) loc j=2
        else loc j=3            
    }    

    /*=========================================================================
                    4: Output
    ===========================================================================*/

    // Case 1: beta
    if "`type'"=="beta" {    
            
        di _n as txt ///
        _col(18) "{hline 4} Treatment Effect Estimate {hline 4}" _n ///
        _col(14) "{c |}" _col(20) "Estimate" _col(39) "Sq. difference" _col(59) "Bias changes" _n ///
        _col(14) "{c |}" _col(37) "from controlled beta" _col(61) "direction" _n ///
        _col(1) "{hline 13}{c +}{hline 64}" _n ///
        _col(1) "Beta" _col(14) "{c |}" _col(18) as result %11.5f `betax' _col(39) %11.3g `distx' _col(63) "`markx'"  _n ///
        _col(1) as txt "Alt. sol. 1" _col(14) "{c |}" _col(18) as result %11.5f `altroot1' _col(39) %11.3g `dist1' _col(63) "`mark1'"  _n ///
        _col(1) as txt "Alt. sol. 2" _col(14) "{c |}" _col(18) as result %11.5f `altroot2' _col(39) %11.3g `dist2' _col(63) "`mark2'"   _n ///
        as txt _col(1) "{hline 13}{c +}{hline 64}"     
    
        di _n as txt ///
        _col(18) "{hline 4} Inputs from Regressions {hline 4}" _n ///
        _col(14) "{c |}" _col(21) "Coeff." _col(39) "R-Squared" _n ///
        _col(1) "{hline 13}{c +}{hline 64}" _n ///
        _col(1) "Constrained" _col(14) "{c |}" _col(18) as res %12.5f `beta_o' _col(39) %5.3f `r_o' _n ///
        _col(1) as txt "Full" _col(14) "{c |}" _col(18) as res %12.5f `beta_tilde' _col(39) %5.3f `r_tilde' _n ///
        _col(1) as txt "{hline 13}{c +}{hline 64}" 

        di _n as txt ///
        _col(18) "{hline 4} Other Inputs {hline 4}" _n ///
        _col(1) "{hline 13}{c +}{hline 64}" _n ///
        _col(1) "R_max" _col(14) "{c |}" _col(18) %5.3f `rmax' _n ///
        _col(1) "Delta" _col(14) "{c |}" _col(18) %5.3f `delta' _n ///
        _col(1) "Models" _col(14) "{c |}" _col(18) "`constrained_model' `full_model'"  _n ///
        _col(1) "{hline 13}{c +}{hline 64}"  

        return scalar beta=`betax'
        return scalar dist=`distx'
        if `j'==2 | `j'==3 {
            return scalar altsol1=`altroot1'
            return scalar altdist1=`dist1'
        }
        if `j'==3 {
            return scalar altsol2=`altroot2'
            return scalar altdist2=`dist2'
        }
        return scalar root_count=`j'
        return scalar delta=`delta'
        
    }

    // Case 2: delta
    if "`type'"=="delta" {
    
        di _n as txt ///
        _col(18) "{hline 4} Bound Estimate {hline 4}" _n ///
        _col(1) "{hline 13}{c +}{hline 64}" _n ///
        _col(1) "delta" _col(14) "{c |}" _col(18) as result %11.5f `boundx' _n ///
        as txt _col(1) "{hline 13}{c +}{hline 64}"     
    
        di _n as txt ///
        _col(18) "{hline 4} Inputs from Regressions {hline 4}" _n ///
        _col(14) "{c |}" _col(21) "Coeff." _col(49) "R-Squared" _n ///
        _col(1) "{hline 13}{c +}{hline 64}" _n ///
        _col(1) "Constrained" _col(14) "{c |}" _col(18) as res %12.5f `beta_o' _col(49) %5.3f `r_o' _n ///
        _col(1) as txt "Full" _col(14) "{c |}" _col(18) as res %12.5f `beta_tilde' _col(49) %5.3f `r_tilde' _n ///
        _col(1) as txt "{hline 13}{c +}{hline 64}" 

        di _n as txt ///
        _col(18) "{hline 4} Other Inputs {hline 4}" _n ///
        _col(1) "{hline 13}{c +}{hline 64}" _n ///
        _col(1) "R_max" _col(14) "{c |}" _col(18) %5.3f `rmax' _n ///
        _col(1) "Beta" _col(14) "{c |}" _col(18) %9.6f `beta' _n ///
        _col(1) "Models" _col(14) "{c |}" _col(18) "`constrained_model' `full_model'"  _n ///
        _col(1) "{hline 13}{c +}{hline 64}"          
            
        return scalar delta=`boundx'
        return scalar beta=`beta'
    }

    return scalar rmax=`rmax'
    return local cmd="reghdfe"
    return local depvar="`depvar_full'"
    return local treatment="`treatment'"
    return local constrained_model="`constrained_model'"
    return local full_model="`full_model'"
    return local type="`type'"
    
    * Restore the full model before exiting
    estimates restore `full_model'
    
end

[Felix Bittmann]

Just FYI, this package might also be interesting for you: https://github.com/mattmasten/regsen...e/vignette.pdf

[Nick Cox]

#11 Thanks for that. It looks as if that code should be in reghdfe_psacalc.ado

help trace

tells you about debugging commands.

[George Ford]

Open the do file and run it entirely. Then try it. Also, here it is saved as ado. I keep this stuff in my ado/personal directory. reghdfe_psacalc.ado

[George Ford]

Here's the demonstration.

xtreg and reghdfe present very different R2, which I suppose explains the small differences. Note that a straight reg with i.id and i.t gives the reghdfe R2, not the xtreg one. I'm sure someone knows why the R2 are so different.

I noticed, however, that psacalc after xtreg uses the within R2, but reghdfe_psacalc is using R2. I'll try to fix that.

clear all

set obs 20
g id = _n
g fe = rchi2(2)
expand 10
bys id: g t = _n

g treated = id <= 5
g control = 1 - treated
g post = t > 5

g x1 = rgamma(5,1)

g y = fe + 0.5*x1 + 2*post*treated + rnormal()

xtset id t
reghdfe y c.post#c.treated, absorb(id t) 
xtreg y c.post#c.treated i.t , fe
xi: xtreg y c.post#c.treated x1 i.t, fe 
psacalc beta c.post#c.treated , mcontrol(_It_2 _It_3 _It_4 _It_5 _It_6 _It_7 _It_8 _It_9 _It_10)

                 ---- Treatment Effect Estimate ----
             |     Estimate           Sq. difference      Bias changes
             |                      from controlled beta    direction
-------------+----------------------------------------------------------------
Beta         |       2.88658                 .285             
Alt. sol. 1  |     -65.09204                 4549             Yes
Alt. sol. 2  |                                                
-------------+----------------------------------------------------------------

                 ---- Inputs from Regressions ----
             |      Coeff.            R-Squared
-------------+----------------------------------------------------------------
Uncontrolled |        1.72188         0.127
Controlled   |        2.35269         0.603
-------------+----------------------------------------------------------------

                 ---- Other Inputs ----
-------------+----------------------------------------------------------------
R_max        |   1.000
Delta        |   1.000
Unr. Controls|   _It_2 _It_3 _It_4 _It_5 _It_6 _It_7 _It_8 _It_9 _It_10
-------------+---------------------------------------------------------------- 

eststo e1: reghdfe y c.post#c.treated , absorb(id t)
eststo e2: reghdfe y c.post#c.treated x1 , absorb(id t)
reghdfe_psacalc beta c.post#c.treated , models(e1 e2)

                 ---- Treatment Effect Estimate ----
             |     Estimate           Sq. difference      Bias changes
             |                      from controlled beta    direction
-------------+----------------------------------------------------------------
Beta         |       2.87917                 .277             
Alt. sol. 1  |     -35.08769                 1402             
Alt. sol. 2  |                                                
-------------+----------------------------------------------------------------

                 ---- Inputs from Regressions ----
             |      Coeff.            R-Squared
-------------+----------------------------------------------------------------
Constrained  |        1.72188         0.522
Full         |        2.35269         0.783
-------------+----------------------------------------------------------------

                 ---- Other Inputs ----
-------------+----------------------------------------------------------------
R_max        |   1.000
Delta        |   1.000
Models       |   e1 e2
-------------+----------------------------------------------------------------