Recently, I’ve spent some time implementing a particular task in both Stata and in Python. Despite this very useful reference, I’ve had a lot of struggle getting the output in both platforms exactly identical. Here is a list of things to watch out for when you’re puzzled by why the two sets of code yield different output.

1. Computing Beta

Suppose you have a dataframe with columns called permno, date, retrf, and mktrf .

• retrf is the excess stock return and mktrf is the excess market return.

The goal is to compute the beta for each stock in a 60-month rolling window with at least 24 valid observations.

Stata Implementation

m = tm(2020m12)
gen byte window = inrange(date,`m'-60,`m'-1)
egen byte obs = count(retrf) if window, by(permno)

egen Mmktrf = mean(mktrf) if retrf<. & window & obs>=24, by(permno)

gen xx = (mktrf-Mmktrf)^2 if retrf<. & window & obs>=24
gen xy = (mktrf-Mmktrf)*retrf if retrf<. & window & obs>=24

egen Mxx = mean(xx), by(permno)
egen Mxy = mean(xy), by(permno)

replace beta = Mxy/Mxx if date ==`m'

Python Implementation

ym = '2020-12-01'
sub_df = df[(df['ldate'] >= ym - pd.DateOffset(months = 60)) &
                 (df['ldate'] <= ym - pd.DateOffset(months = 1))]

covariance = sub_df[['ldate', 'daret_rf', 'mktrf', 'permno']].groupby('permno').cov(min_periods = 24, ddof = 0).reset_index()
numerator = covariance[['permno', 'daret_rf']][1::2].set_index('permno')['daret_rf']
denominator = covariance[['permno', 'mktrf']][1::2].set_index('permno')['mktrf']

betas =  DataFrame(numerator.divide(denominator)).reset_index()

Nota Bene

• The above Stata implementation uses population covariance (as seen above) while the Python covariance uses the sample covariance.
• The above Stata implementation drops mktrf if retrf is missing, but Python does not do this.
  • It turns out that this discrepancy leads to sensitive estimates of beta, especially when the number of observations is close to the minimum required threshold.

2. Generating Cuts

Suppose you have a dataframe with columns called permno, signal , and date and you want to create 10 equally sized portfolios at each date based on the signal value.

Stata Implementation

egen portfolio = xtile(signal), by(date) nq(10)

Python Implementation

df.groupby('date')['signal'].transform(lambda x : pd.qcut(x, 10))

Nota Bene

• There’s a slight difference in how Python and Stata assign the portfolios.
• For example, consider a stock that ranks 294th out of 997 stocks at date = 1990.07.
  • Note that 0.3 * 97.7 = 293.1
• In Stata, this stock is classified as portfolio #3, while in Python it’s classified as portfolio #4.

3. Winsorizing Variables

Consider a dataframe with columns called beta and date and winsorizing beta at the 2.5% level at each date.

Stata Implementation

winsor2 beta, replace cuts(2.5 97.5) by(date)

Python Implementation

df.groupby('date')['beta'].transform(lambda x : x.clip(lower = x.quantile(0.025),
                                                                                  upper = x.quantile(0.975))

4. Replacing Multiple Observations with Median Values

Consider a dataframe with columns permno, signal, and date and suppose there are multiple observations with the same permno-date pair in this dataset. We are then interested in replacing these duplicates with their median values.

Stata Implementation

egen Msignal = median(signal), by(permno date)
replace signal = Msignal

Python Implementation: This is somewhat more convoluted —

duplicates = df['signal'].groupby([df['permno'], df['date']]).count() \
                  + df['signal'].isnull().groupby([df['permno'], df['date']]).sum().astype(int)
medians = df[['permno', 'date', var]].groupby(['permno', 'date']).median()['signal']

for ind in duplicates[duplicates > 1].index:
   sub_df = df['signal'][(df_comp_q['permno'] == ind[0]) & (df['date'] == ind[1])]
  if len(sub_df.dropna()) > 0:
     median_value = np.nanmedian(sub_df)
    df[var][(df['permno'] == ind[0]) & (df['date'] == ind[1])] = median_value

5. PCA

You can do a PCA on the correlation matrix or the covariance matrix. Suppose you have 20 columns of data called col1 ~ col20 in a dataframe called df . You’re interested in the first 4 PCs.

Stata Implementation

pca col*, corr
predict pca1 pca2 pca3 pca4, score

Python Implementation: Using PCA Module

from sklearn.decomposition import PCA
from sklearn.preprocessing import StandardScaler
pca = PCA(n_components = 4)
sc = StandardScaler()

components = DataFrame(data = pca.fit_transform(sc.fit_transform(df)), columns = ['PC1', 'PC2', 'PC3', 'PC4'], index = df.index)

# Sanity Check
components.corr()

• This implementation applies the obtained PC coefficients to the scaled matrix, not the original matrix. (In the newest version, this is what Stata seems to do).

Python Implementation: Using PCA Module

corr_matrix = df.corr()

components = {}
eig_vals, eig_vecs = np.linalg.eig(corr_matrix)
components ['PC1'] = df.dot(eig_vecs[:,0])
components ['PC2'] = df.dot(eig_vecs[:,1])
components ['PC3'] = df.dot(eig_vecs[:,2])
components ['PC4'] = df.dot(eig_vecs[:,3])
components = DataFrame(components )

• This implementation computes the PC from the correlation matrix explicitly and then applies the obtained coefficients to the original matrix.

6. Using Stata in Python

With Stata17, you can now use Python in Stata. Here’s the preamble:

import stata_setup
stata_setup.config('C:\Program Files\Stata17', 'se')
import sys
sys.path.append('STATA_SYSDIR/utilities')
from pystata import stata

The link here contains example usages.

Example Usage

stata.pdataframe_to_data(df, force = True)
stata.run('summarize')

7. Running Regressions

While there are many ways to run OLS regressions in Python, I find it easiest to use R-style formulas. For more information, please check this patsy documentation.

• Simple OLS Regression

  m1 = smf.ols(formula = 'y ~ x', data = df).fit(use_t = True)
  

• Heteroscedasticitiy-consistent Standard Errors

  m1 = smf.ols(formula = 'y ~ x', data = df).fit(cov_type = 'HC1', use_t = True)
  

• Clustering Standard Errors by One Dimension Note: Make sure to drop missing values since clustering may result in misaligned rows.

  m1 = smf.ols(formula = 'y ~ x', data = df).fit(cov_type = 'cluster',
                                                                            cov_kwds = {'groups' : df['firm'],
                                                                            use_t = True)
  

• Clustering Standard Errors by Two Dimensions Note: Similar as before, but one needs to encapsulate it in numpy.array()

  m1 = smf.ols(formula = 'y ~ x', data = df).fit(cov_type = 'cluster',
                                                                            cov_kwds = {'groups' : np.array(df[['firmid', 'year']]),
                                                                            use_t = True)
  

• Newey-West Adjustment for Standard Errors Note: The use_correction : True is necessary to equate it with output from Stata.

  m1 = smf.ols(formula = 'y ~ x', data = df).fit(cov_type = 'HAC',
                                                                            cov_kwds = {'maxlags' : 3,
                                                                                              'use_correction' : True},
                                                                            use_t = True)
  

Other Resources:

• Diff-in-diff estimators in Python (link to Jupyter notebook)

8. Constructing Weighted Averages

Suppose you want to construct weighted averages of r_it across i for a given t using a variable called wealth. The following code generates such weighted averages.

df['r_St'] = df.groupby(['t']).apply(lambda x : np.average(x['r_it'], 
                                                                                          weights = x['wealth'])

Setting weights equal to None constructs an equal-weighted average.

9. Binscatter

The easiest way to implement a binscatter in Python is to use the seaborn library.

sns.regplot(data = df_plot, x = x_var, y = y_var, 
                  x_bins = 100, ci = None, fit_reg = False, color = 'gray')

10. Balancing Panel

The easiest way to balance a panel in Python is using the unstack and the stack comment:

df = df.set_index(['firmID', 'year'].unstack().stack(dropna = False)