Methods for measurement of delays in the bankruptcy process

by Dhananjay Ghei and Shubho Roy.

When dealing with distress, the single big idea for avoiding value destruction is speed of action. Delays destroy value. When we cope with failing firms (through a bankruptcy code) or failing banks (through a resolution corporation), the measure of institutional quality is the recovery rate, and the key determinant of the recovery rate is speed of action.

Methods for measurement of delays

The conventional method for assessing delays is to run a survey of insolvency practitioners (judges, lawyers, accountants). In this article, we offer a new idea for measuring one part of the delay.

There are many parts of the bankruptcy process. Many societies try to avoid the problem of failed firms by putting things off. This leads to a delay between the actual firm failure and the commencement of the formal bankruptcy process. Could this be measured directly? We could start from the date of the end (the date of the official legal action about firm closure), and look back to the date on which extreme credit distress is identified. The gap between the two events will serve as a proxy for delays of the bankruptcy process.

Let's play this idea out with the failure of Kingfisher Airlines.

The decisive end date is 18th November, 2016, when the Karnataka High Court ordered the winding up of Kingfisher Airlines (See here).

Conventional notion of delay: The case was filed by an unsecured (foreign) creditor on 18th September, 2012. It took the legal system four years to come to the conclusion that Kingfisher Airlines is insolvent.

Looking back using interest coverage ratio (ICR) as a measure of credit stress. This is defined as the ratio of the earnings of a company (EBITDA) and the total interest payments due. The logic being that if a company's income is not enough to even pay the interest on its borrowing, there is acute distress. The literature in finance uses ICR in a couple of different ways when assessing firm distress. Here are three examples:

ICR below 1

Distress is where the ICR falls below 1 in any year. Example: Claessens, Djankov, and Lang 1998..

ICR below 0.75

Distress is where the ICR for any financial year falls below 0.75. Example: Love, 2010.

P-ICR

Distress is Persistent-ICR, which is defined as an ICR of less than 1 for three consecutive financial years. Example: Chung and Ratnovski, 2016

We extract the annual financial data for Kingfisher Airlines from CMIE Prowess database. The data is available from March 2001 till March 2013. Figure 1 below shows when Kingfisher Airlines met these criteria. It broke into ICR-below-1 and ICR-below-0.75 in March 2005 (red line) and it achieved P-ICR in March 2007 (orange line). The court finally ordered winding up 3520 days after Kingfisher meeting all the three criteria (gray line). The black line (October 20, 2012) is when DGCA suspended Kingfisher's license.

Figure 1: Interest cover ratio for Kingfisher Airlines

Looking back using Distance to Default as a measure of credit stress. Distance to Default (DtD) is measured as the difference between the asset value of the firm and the face value of its debt, scaled by the standard deviation of firm's asset value. It measures the distance (in standard deviations) between the expected value of the firm and the "default point" (face value of the debt). Thus, lower values of DtD imply that the firm is more likely to default on its financial obligations. DtD is implemented in R using the ifrogs package developed by IGIDR Finance Research Group.

Figure 2: Distance to default for Kingfisher Airlines

Figure 2 shows the DtD for Kingfisher Airlines from June 2007 onwards. As the listing took place on 12 June 2006, the DtD calculation only commences on 12 June 2007. The red line shows the date when DtD was below 1 for the first time (August 1, 2008). DtD reached its lowest value of 0.30 on July 6, 2009 (orange line). Thus, at that point the likelihood of default was the highest (19 percent). 

Conclusions

The approach shown above with Kingfisher is potentially scalable into large datasets and all countries. This could yield large-scale measurement, at the level of individual bankruptcy transactions, about delays. This can yield useful summary statistics.

The new insolvency law should reduce the time between financial insolvency and the legal recognition of the same. This should yield measurable gains when these techniques are applied to the data in the future.

Our approach yields micro data about delays which can be analysed in order to explore cross-sectional and time-series variation. Such analysis is not feasible with the conventional measurement of delays that is based on surveying practitioners.



Dhananjay Ghei and Shubho Roy are researchers at the National Institute for Public Finance and Policy. The authors would like to thank Radhika Pandey, Rajat Kochhar and Mohit Desai for their inputs.

Author: Dhananjay Ghei

Dhananjay Ghei is a researcher at the National Institute for Public Finance and Policy. On this blog:

• Methods for measurement of delays in the bankruptcy process, 25 November 2016. 
• The Diwali effect in Delhi air quality, 28 October 2016.
  
• Describing Delhi's air quality crisis, 13 October 2016.
  
• The gains to US GDP from a Doing Business score of 100, 11 July 2016.
• Sophisticated clustered standard errors using recent R tools, 15 June 2016.
  

The Diwali effect in Delhi air quality

by Dhananjay Ghei, Arjun Gupta and Renuka Sane

As Diwali approaches, we have learned to worry about air quality. Over the last few years, several studies have noted the increase in pollution levels during the period of Diwali owing to increase in commercial activity and firework displays. However, as we show in our previous article, there is considerable variation in PM 2.5 levels in Delhi in terms of location/time/month:

1. Time Effect: The effect of diwali is not uniform throughout the day and is more prevelant at particular time of the day than other times. We also need to adjust for the confounding effect of time: pollution levels are high during the night and low during the day.
2. Location Effect: Several areas of Delhi are severly polluted throughout the time, whereas others see large variations in their pollution levels. All these reasons make it difficult to attribute the entire increase in PM2.5 on Diwali.
3. Month Effect: The day of Diwali Festival varies in the Gregorian Calendar between the 17th October and 15th November every year. Existing pollution levels are already high when compared to the annual average. This is a confounding effect.

It is possible that the bad air that we see in Delhi at the time of Diwali is just the bad air quality in winter, and is not causally impacted upon by Diwali. In this article, we attempt to quantify the increase in the PM 2.5 levels during the Diwali period. Does Diwali have an impact upon air quality? If so, by how much?

Issues in research design

The opportunity to identify a Diwali effect comes from the fact that Diwali is a `moving holiday' which takes place on a different day of each year. If this were not the case, it would be strongly correlated with changing climate.

Our ability to analyse these questions is greatly hampered by the lack of data. As of today, the data only runs from 1/2013 to 10/2016.

The air pollution caused by fireworks includes many contaminants. The data that we are studying covers only pm2.5.

Pollution levels on Diwali

The data used for the analysis comes from the US Consulate based in Chanakyapuri and the Central Pollution Control Board for 4 locations (R K Puram, Punjabi Bagh, Mandir Marg, Anand Vihar). The data consists of hourly PM 2.5 levels across the five locations from January 2013 to October 2016. We winsorise the data at 1% on both ends to remove the extreme tail values.

The effect of Diwali on pollution levels

We first estimate the effect of Diwali on daily data using an event study. We aggregate the hourly concentration of PM2.5, at each location, to arrive at the daily numbers. The day of the Lakshmi Puja is taken as the event day. Therefore, we get 3 events for each location. Next, we calculate the percentage change in PM2.5 concentration levels by differencing the logarithm of PM2.5 values. These are then re-indexed to show the cumulative change over a 20 day window.

Event study showing the change in PM2.5 around Diwali date (in days)

The solid line represents the average cumulative percentage change in PM2.5 values during the window, whereas the dashed line represents the confidence intervals calculated using the bootstrapped standard errors. We see that pollution levels start increasing one day before Diwali, and increase till two days after Diwali. It is also interesting to note that the increase in the pollution levels is significant during the two days after Diwali. This can be attributed to the fact that Diwali celebrations begin only on the night of Diwali, thereby leading to a significant increase the next day, as well as Diwali being celebrated over an extended period of time.

We now come at the same set of questions using a regression.

Contribution of Diwali on PM2.5: Regression analysis

Since Diwali is celebrated over a number of days we also define the following models:

1. Diwali=t: Diwali
2. Diwali={t-1:t+1}: 3 Days (day before Diwali, Diwali, day after Diwali)
3. Diwali={t-1:t+2}: 4 Days (preceding day to two days after Diwali)

The model is as follows:

\[ PM2.5_{it} = \alpha + \beta_1*Diwali_{t}+ \beta_2*Diwali_{t}*l_{i} + m_t + h_t + l_i+\epsilon_{it} \]

where, $i$ is location, and $t$ is time. Here, PM 2.5 is the hourly measured levels of the pollutant. The first model takes Diwali to be only the date of Diwali, second model defines the Diwali days from one day before to one day after and the third model considers Diwali from the preceding day to two days after Diwali. In addition, we have month ($m_t$), location ($l_i$), and hour ($h_t$) fixed effects. The base for the location interaction term is Anand Vihar. Robust standard errors are used for our analysis throughout.

	Dependent variable:
	Hourly PM2.5 Concentration
	Diwali=t	Diwali={t-1:t+1}	Diwali={t-1:t+2}
	(1)	(2)	(3)
Diwali	-3.720	98.687	134.709
	t = -0.177	t = 8.496***	t = 13.181***
Chanakyapuri*Diwali	17.270	-75.878	-87.035
	t = 0.638	t = -5.100***	t = -6.692***
Mandir Marg*Diwali	73.078	-67.943	-66.844
	t = 2.606***	t = -4.450***	t = -4.979***
Punjabi Bagh*Diwali	65.630	-49.033	-52.254
	t = 2.374**	t = -3.254***	t = -3.945***
R K Puram*Diwali	63.348	-54.228	-67.094
	t = 2.291**	t = -3.589***	t = -5.055***
Month FE	Yes	Yes	Yes
Location FE	Yes	Yes	Yes
Hour FE	Yes	Yes	Yes
Observations	118,847	118,847	118,847
R2	0.264	0.264	0.266
Adjusted R2	0.264	0.264	0.266
F Statistic (df = 39; 118803)	1,091.020***	1,094.274***	1,103.673***

The first model (Column 1) shows that the baseline effect (i.e. at Anand Vihar) is not statistically different from non-Diwali days. For locations, other than Chanakyapuri, there is a differential effect on Diwali relative to Anand Vihar on Diwali. For instance, Diwali adds on an average 69.35 (73.07-3.72) µg/m³ PM2.5 particulate matter in air at Mandir Marg relative to Anand Vihar.

When we consider the second (Column 2) and third (Column 3) specifications, there is a statistically significant effect in Anand Vihar. The average particulate matter is 99 µg/m³ higher when we consider a two day Diwali, and 135 µg/m³ when we consider a three day Diwali period. While this may not seem much, given the already degraded air quality during these months, Diwali makes the pollution level reach alarming levels (>400, the monthly average in October November is around 340) which can have severe impacts on the health of people.

The Diwali effect is lower in other other locations relative to Anand Vihar. Thus, we see, that on the main day of Diwali, Anand Vihar is not too different from other days, while other locations have more pollutants relative to Anand Vihar. However, once we take into account 1-2 days after Diwali, we see that Anand Vihar is the most polluted location, and other locations have lower pollutants relative to Anand Vihar.

Conclusion

Very little is known, at present, about air quality and Diwali. Using the admittedly weak data resources, we have begun analysing this question here.

To the extent that these results are persuasive, they could help individuals plan strategies to avoid being in Delhi on these days. There is also a case for a Pigouvian tax on fireworks, in order to overcome the externality.

Previous work on Diwali, which helps us see other dimensions of Diwali, includes: Seasonal adjustment with Indian data: how big are the gains and how to do it by Rudrani Bhattacharya, Radhika Pandey, Ila Patnaik, Ajay Shah, and IEDs in Diwali and Toxic chemicals in Holi by Ajay Shah.

Describing Delhi's air quality crisis

by Dhananjay Ghei, Arjun Gupta and Renuka Sane.

One of the most important elements of public health is regulatory interventions that yield clean air.  In late 2016, we await the air quality crisis of the Delhi winter with trepidation. A few attempts at solving the problem have begun. The Government of Delhi experimented with an odd-even policy to regulate traffic between 1 January 2016 to 15 January 2016, and then between 15 April 2016 to 22 April 2016. The results of these experiments have been mixed [here and here].

What you measure is what you can manage. Only when we are able to marshal evidence in a systematic way about the extent and nature of the problem, will we be able to design and deliver a response. The measurement of air pollution in Delhi has begun on a small scale. In this post, we describe patterns seen in the available data.

Why is PM 2.5 a good measure?

There are many pollutants in the air such as carbon monoxide (CO), nitric oxide (NO), nitrogen dioxide (NO2), ozone (O3). The worst among these is small particulate matter, or PM 2.5, which are a mixture of solid and liquid droplets floating in the air whose diameters are less than 2.5 micrometers. These fine particles are produced from all types of combustion, including motor vehicles and power plants and some industrial processes.

The health impact from pollution is a complex transform of exposure to all pollutants. However, of the pollutants, PM 2.5 particles are considered the most harmful as they are able to enter deep into the respiratory tract, reaching the lungs. This can cause short-term health effects such as eye, nose, throat and lung irritation, coughing, sneezing, runny nose and shortness of breath, and in the long-term can affect lung function and worsen medical conditions such as asthma and heart disease. We, therefore, narrow our attention to the measure of PM 2.5. The unit of measurement of PM 2.5 is µg/m³ and the breakpoints of raw PM 2.5 values by the US Environmental Protection Agency are the following:

24-hr PM 2.5	AQI Categories	Health Effects Statements
0.0-12.0	Good	None
12.1-35.4	Moderate	Respiratory symptoms possible in unusually sensitive individuals, possible aggravation of heart or lung disease in people with cardiopulmonary and older adults.
35.5-55.4	Unhealthy for Sensitive Groups	Increasing likelihood of respiratory symptoms in sensitive individuals, aggravation of heart or lung disease and premature mortality in people with cardiopulmonary disease and older adults.
55.5-150.4	Unhealthy	Increased aggravation of heart or lung disease and premature mortality in people with cardiopulmonary disease and older adults; increased respiratory effects in general population.
150.5-250.4	Very Unhealthy	Significant aggravation of heart or lung disease and premature mortality in people with cardiopulmonary disease and older adults; significant increase in respiratory effects in general population
250.5-500	Hazardous	Serious aggravation of heart or lung disease and prematuremortality in people with cardiopulmonary disease and older adults; serious risk of respiratory effects in general population.

 

Data

We fetch raw PM 2.5 values from two data sources on pollution in Delhi. The first is put out by the US Embassy based in Chanakyapuri. In addition, the Central Pollution Control Board also puts out real time data for various locations across India. We select 4 locations which provided us with the most consistent dataset. This gives us a total of 5 locations for which we have data:

1. R K Puram
2. Punjabi Bagh
3. Mandir Marg
4. US Embassy (Chanakyapuri)
5. Anand Vihar

We use hourly data from the locations mentioned above for a time period from January 2013 to October 2016. It should be noted that values are missing from certain sections of the data. These missing observations are excluded from our analysis.

Drawing upon the Chinese experience, it's interesting to ask: Do the Indian government sources tally with the US Embassy data?  We can't say, as there is no measurement for a location near the US Embassy by the CPCB.

Dimensions of variation

These are three types variations seen in PM 2.5.

Figure 1: Variation by time of day

Time Effect: Figure 1 above shows the variation in hourly pollution levels during different days of a week. Darker colors represent increased PM 2.5 matter in the air. We see that the pollution levels are low during the day, but start increasing post 6 p.m. and remain elevated till 9 a.m. of the next day. The average PM 2.5 concentration from 6 p.m. to 9 a.m. is 140 µg/m³, whereas the average PM 2.5 concentration from 9 a.m. to 6p.m. is 108 µg/m³. PM 2.5 levels in the range of 101-200 can cause breathing discomfort to anyone with prolonged exposure to the air during these times. This graph suggests that a measure that restricts traffic during the day such as the odd-even policy is unlikely to be as effective as a measure that restricts emissions at night.

Figure 2: Variation by month

Month Effect : Figure 2 shows the hourly variation in pollution levels during different months of the year. Note that the scale for this figure is different from that used in Figure 1. The monsoon months have the lowest levels of PM 2.5 particulate matter. Larger particles are settled in few hours due to gravity, but smaller particles such as PM 2.5 are removed by precipitation. Winters have the highest levels of PM2.5 matter in the air, on account of low wind speed and high relative humidity. PM 2.5 concentration reaches above 200 in the winter months, which can cause respiratory illness to people on prolonged exposure and puts people with respiratory illness, and heart disease on a far greater risk.

Figure 3: Variation by location

Location Effect: Figure 3 shows the hourly variation in pollution levels at the five locations where instruments are available. Chanakyapuri seems to perform better than other areas of Delhi, in terms of PM 2.5 particulate matter. Anand Vihar has the highest pollution levels amongst the 5 different locations, and has severe levels of air pollution in the night. This can cause respiratory impact even on healthy people, and serious health impacts on people with lung/heart diseases.

Thus, we see that there is a strong location effect on pollution levels. This can be due to the varying population densities of these locations as well as the proximity to industries etc. This could lead to location-specific policy initiatives such as closing down factories or modifying vehicular traffic.

Reproducible research

Data and R code.



Dhananjay Ghei and Arjun Gupta are researchers at the National Institute of Public Finance and Policy. Renuka Sane is an academic at the Indian Statistical Institute, Delhi Centre.

The gains to US GDP from a Doing Business score of 100

by Dhananjay Ghei and Nikita Singh.

Can a country achieve growth by implementing large pro-business reforms? If yes, then how much growth is really possible from such reforms? In a recent WSJ op-ed, Cochrane takes a stab at this question for the United States. Using data from the World Bank's ease of doing business index, Cochrane claims there is a log-linear relationship between GDP per person and business climate. By extrapolating this relationship out of sample, he predicts that the US would register a 209% improvement in per capita income (or, 6% additional annual growth if the required reforms are implemented over the next 20 years) by achieving the ease of doing business index value of 100.

Brad Delong disagrees. He fits a fourth-degree polynomial on the same data. He justifies this on the grounds that the third degree coefficient is negative and statistically significant. His forecast shows that an increase in the index value beyond 90 would actually lead to a lower GDP per person. Figure 1 juxtaposes the log-linear and polynomial regression fit, and we can see how the two views are sharply different. The straight line yields higher and higher GDP as you go to 100; the polynomial droops off at the end.

Figure 1: The analysis of John Cochrane and Brad DeLong

Areas of concern

There are many areas of concern with this analysis:

1. Assuming linearity is surely a stretch. But polynomial regressions are a bad way to deal with nonlinearity. In particular, polynomial regressions are very fragile at the end points. This can be easily seen in Figure 2 as the prediction interval increases at edges of the data. In addition, extrapolation using a polynomial is almost always sure to give a wrong answer as the curvature of the polynomial is unidentified outside the sample.
2. Using a cross-sectional regression with one variable is a poor guide to the causal relationships. Labour and capital matter to GDP per capita. There are stark differences in law and governance, institutions and culture across countries; it is unlikely that the doing business score is a sufficient statistic.
3. Hallward-Driemeier and Pritchett (2015) show that the "doing business index" is not a good reflection of how the laws on paper are implemented in reality. The main point of their argument is that better de jure regulations do not necessarily imply improved de facto outcomes specially when a country has weak governmental capabilities for implementation and enforcement. Even if the US does well on the rule of law, and this gap between rules and deals is absent, this is a serious issue for many (most?) observations in the dataset.

Figure 2: The 95% prediction interval for the polynomial regression

Can we do better?

Criticisms 2 and 3 are hard to handle. But a little bit of statistics helps us do better on the first. We use non-parametric regression as a way to have nonlinearity in the relationship between business climate and GDP per person without having to take a stand on a particular functional form. This involves three steps:

1. Selecting an optimum bandwidth using cross-validation
2. Estimating a nonparametric model using the chosen bandwidth
3. Tests of statistical significance and specification

We use a second order Gaussian kernel and fit a local linear estimator to identify the functional form in sample. Business climate is significant at 1% level in the local linear non-parametric model. Moreover, based on a lower cross validation score, the non-parametric regression is favoured.  In addition, we do a bunch of robustness tests by changing the type of kernel and regression. The results do not change much in either of the cases. These calculations were done in  R using the np package.

Figure 3: Non-parametric regression gives us the best of both worlds

The results, shown above, show that there is nonlinearity in the data. The linear model used by Cochrane is not appropriate. But we're better off as compared with using a polynomial regression; the confidence interval is tighter at the edges.

Figure 4 superposes the three models. The coloured dots show the predicted value of GDP per person using the three different specifications when the doing business index takes the value of 100.

Figure 4: Comparing the three predictions

Our nonparametric estimate shows that gains from achieving a score beyond 90 are increasing and somewhere in between Cochrane and DeLong's numbers. Cochrane predicts that the US would achieve 6% additional annual growth for 20 years by moving to a score of 100. If we go out of sample to estimate using the nonparametric fit, this shows an annual growth of 2.22% for the next 20 years. This is not something to laugh at, but it's a smaller, and we think a more plausible estimate.

References

Hallward-Driemeier, Mary and Lant Pritchett. 2015. "How Business Is Done in the Developing World: Deals versus Rules." Journal of Economic Perspectives, 29(3): 121-40.

Tristen Hayfield and Jeffrey S. Racine (2008). Nonparametric Econometrics: The np Package. Journal of Statistical Software 27(5). URL http://www.jstatsoft.org/v27/i05/.


Dhananjay Ghei is a researcher at the National Institute of Public Finance and Policy. Nikita Singh is a MRes. student at London School of Economics and Political Science. The authors thank Ajay Shah for valuable discussions and feedback.

Sophisticated clustered standard errors using recent R tools

by Dhananjay Ghei

Many blog articles have demonstrated clustered standard errors, in R, either by writing a function or manually adjusting the degrees of freedom or both (example, example, example and example). These methods give close approximations to the standard Stata results, but they do not do the small sample correction as the Stata does.

In recent months, elegant solutions have come about in R, which push the envelope on functionality, and yield substantial improvements in speed. I use the test dataset of Petersen which is the workhorse of this field.

The problem

In regression analysis, getting accurate standard errors is as crucial as obtaining unbiased and consistent estimates of the regression coefficients. Standard errors are important in determining the accuracy of the coefficients and thereby, affecting hypothesis testing procedures.

The correct nature of standard errors depends on the underlying structure of the data. For our purposes, we consider cases where the error terms of the model are independent across groups but correlated within groups. For instance, studies with cross-sectional data on individuals with clustering on village/state/hospital level. Another example could be difference in difference regressions with clustering at a group level. Clustered standard errors allow for a general structure of the variance covariance matrix by allowing errors to be correlated within clusters but not across clusters. In such cases, obtaining standard errors without clustering can lead to misleadingly small standard errors, narrow confidence intervals and small p-values.

Clustered standard errors can be obtained in two steps. Firstly, estimate the regression model without any clustering and subsequently, obtain clustered errors by using the residuals. Clustered standard errors can be estimated consistently provided the number of clusters goes to infinity. However, the variance covariance matrix is downward-biased when dealing with a finite number of clusters. One of the methods commonly used for correcting the bias, is adjusting for the degrees of freedom in finite clusters.

R and Stata codes

The code below shows how to compute clustered standard errors in R, using the plm and lmtest packages. Petersen's dataset can be loaded directly from the multiwayvcov package. Pooled OLS and fixed effect (FE) models are estimated using the plm package.

# Loading the required libraries
library(plm)
library(lmtest)
library(multiwayvcov)

# Loading Petersen's dataset
data(petersen)
# Pooled OLS model
pooled.ols <- plm(formula=y~x, data=petersen, model="pooling", index=c("firmid", "year")) 
# Fixed effects model
fe.firm <- plm(formula=y~x, data=petersen, model="within", index=c("firmid", "year")) 

Clustered standard errors can be computed in R, using the vcovHC() function from plm package. vcovHC.plm() estimates the robust covariance matrix for panel data models. The function serves as an argument to other functions such as coeftest(), waldtest() and other methods in the lmtest package. Clustering is achieved by the cluster argument, that allows clustering on either group or time. The type argument allows estimating standard errors by allowing for heteroskedasticity across groups. Recently, the plm package introduced the small sample correction as an option to the "type" argument of vcovHC.plm() function. This is switched on by specifying type="sss".

# OLS with SE clustered by firm (Petersen's Table 3)
coeftest(pooled.ols, vcov=vcovHC(pooled.ols, type="sss", cluster="group"))  

# OLS with SE clustered by time (Petersen's Table 4)
coeftest(pooled.ols, vcov=vcovHC(pooled.ols, type="sss", cluster="time")) 


# FE regression with SE clustered by firm
coeftest(fe.firm, vcov=vcovHC(fe.firm, type="sss", cluster="group")) 

# FE regression with SE clustered by time
coeftest(fe.firm, vcov=vcovHC(fe.firm, type="sss", cluster="time")) 

Stata makes it easy to cluster, by adding the cluster option at the end of any routine regression command (such as reg or xtreg). The code below shows how to cluster in OLS and fixed effect models:

webuse set http://www.kellogg.northwestern.edu/faculty/petersen/htm/papers/se/
webuse test_data.dta, clear

* OLS with SE clustered by firm (Petersen's Table 3)
reg y x, vce(cluster firmid)
* OLS with SE clustered by time (Petersen's Table 4)
reg y x, vce(cluster year)

* Declaring dataset to be a panel
xtset firmid year
* FE regression with SE clustered by firm
xtreg y x, fe vce(cluster firmid)
* FE regression with SE clustered by time
xtreg y x, fe vce(cluster year) nonest

The table given below shows a comparison of the standard errors computed by R and Stata. The standard errors computed from R and Stata agree up to the fifth decimal place.

Model	SE (in R)	SE (in Stata)
OLS with SE clustered by firm	0.05059	0.05059
OLS with SE clustered by time	0.03338	0.03338
FE regression with SE clustered by firm	0.03014	0.03014
FE regression with SE clustered by time	0.02668	0.02668

Performance comparison

I run benchmarks for comparing the speed of Stata MP and R for each of these models on a quad-core processor. The results show that R is faster than Stata. In order to do parallelisation, I set the number of processors that Stata MP will use as 4. An example of the benchmarking code in Stata is given below:

* Stata benchmarking program : Example
set processors 4
timer clear
timer on 1
bs, nodrop reps(1000) seed(1): reg y x
timer off 1
timer list

Parallelisation in R is done using standard R packages. An example of the benchmarking code in R is given below:

# R benchmarking program : Example
library(doParallel)
library(rbenchmark)
set.seed(1)
c <- detectCores()
cl <- makeCluster(c)
ols.benchmark <- mcparallel(benchmark(lm(y~x, petersen), replications=1000))
mccollect(ols.benchmark)
stopCluster(cl)

The table below shows a comparison of R and Stata MP for each of these models. The average time is calculated as the ratio of elapsed time to the number of replications. Relative efficiency is defined as the ratio of the average time taken by Stata MP to the average time taken by R. It turns out that the R is faster.

Model	Replications	Average time (R - 4 core)	Average time (Stata MP - 4 core)	Relative efficiency
OLS with SE clustered by firm	1000	0.0737	0.1635	2.22
OLS with SE clustered by time	1000	0.0557	0.0742	1.33
FE regression with SE clustered by firm	1000	0.0880	0.3176	3.61
FE regression with SE clustered by time	1000	0.0729	0.1118	1.53

Multi-level clustering in R

Two way clustering does not have a routine estimation procedure with most of the Stata commands (except for ivreg2 and xtivreg2). There are a few codes available online (See for example, here and here) that do two way clustering. This is easily handled in R, using the vcovDC.plm() function. The function can be used in a similar fashion as vcovHC.plm().

# OLS with SE clustered by firm and time (Petersen's Table 5)
coeftest(pooled.ols, vcov=vcovDC(pooled.ols, type="sss"))

A more recent addition, multiwayvcov package is useful for clustering on multiple levels and, in computing bootstrapped clustered standard errors. The package supports parallelisation thereby, making it easier to work with large datasets. Two functions are exported from the package, cluster.vcov() and cluster.boot(). cluster.vcov() computes clustered standard errors, whereas, cluster.boot() calculates bootstrapped clustered standard errors. The code for replicating Petersen's results is available in the reference manual of the package. One limitation of cluster.vcov() is its inability to work with plm objects. This is because the package imports estfun() from the sandwich package, which is not compatible with plm objects.

R code

Here's the R code to reproduce the results.

Dhananjay Ghei is a researcher at the National Institute of Public Finance and Policy. He thanks Ajay Shah, Vimal Balasubramaniam and Apoorva Gupta for valuable discussions and feedback.