Clustering SE at 2 levels I saw some papers that cluster standard errors at two levels (let's say firm and year) and I found R and STATA packages that can do it. But what is the idea behind clustering at two levels? And how can I state that observations between clusters are still iid this way?
 A: The easiest way to understand this is with a compound error term. To make things a bit more concrete, let's say $y$ is the income of teachers in various schools who specialize in teaching different subjects:
$$y_{ish} = x_{ish}'\beta + \alpha_{is} + \delta_{ih} + \varepsilon_{ish},$$
where $i$ indexes teachers, $s$ indexes schools, and $h$ indexes subjects (math, biology, physics, science, gym, basket-weaving). Observable characteristics like education, field and tenure that vary across teachers are in $x$. The last three terms contain variables unobserved in the data, where the last one is just an idiosyncratic error with the standard assumptions.
The two important properties of the errors are their variance and the covariance/correlation within and across clusters. With $n$ observations, you have an $n \times n$ variance-covariance matrix that looks like this:

The variance $\sigma^2$s along the diagonal are all the same (homoskedastic), and the covariances (off-diagonal terms) are all zero. This says that an error that makes teacher A's wage higher is uncorrelated with the error for teacher B, and so on for all pairs. The variance of the errors is also the same for all teachers. This is what the variance-covariance matrix for the $\varepsilon$ part looks like.
Clustered SEs allow some off-diagonal terms to be non-zero as long as observations belong to the same cluster. They also allow the variances to be different. This doesn't impose a lot of other restrictions on what that might look like.
What are some examples?
The $\alpha_i$s are allowed to be correlated within schools but not across them. For example, some schools have principals with a STEM background, and they may reward STEM teachers with more money. That means the STEM teachers in such schools will all have positive salary errors, and the non-STEM teachers will have correlated negative errors. But the STEM principals cannot alter salaries at different schools, so you don't expect the $\alpha_i$s to be correlated across schools, only within them. You can also have a multitude of other factors operate like this here and at other schools, so it's not just a matter of adding school fixed effects in your model. Moreover, you might also have some unionized schools, where the variance of the errors is smaller because the union bargained for more rigid wage policies that limit inequality. Other schools may have performance-based pay so that the errors have higher variance and less within-school correlation, given some teacher characteristics. In short, all these things can make the variance-covariance matrix more complicated.
The $\delta_i$s are field-specific. For example, let's say that there is a perceived shortage of art teachers, and there is a national program that subsidizes salaries for them, indexed to a local cost of living index. That will create positive errors that are correlated for all teachers whose field is artsy. Or you might expect there to be more variation in some fields than others. A football coach at a high school in TX where the whole town comes out to watch the game on Friday night may have a very high salary, but there just isn't an equivalent of that for math, so the variance in math salaries will be much smaller than for gym teachers.
The variance-covariance matrix with one-way clustering looks more like this:

The $\Sigma$s sub-matrices may not be identical and may have some arbitrary off-diagonal elements. Here's what a four-teacher school submatrix might look like:

In the big version, most of the off-diagonal elements would be zero to reflect that there is no across-school correlation. It will look like a chessboard with uneven shading on the dark blocks and an irregular grid.
With k-way clustering, you would have $k+1$ such matrices. You could use the formula for the variance of random numbers to collapse to just one matrix by calculating the $Cov(\alpha_{is} + \delta_{ih} + \varepsilon_{ish},\alpha_{js} + \delta_{jh} + \varepsilon_{jsh})$, but that is not the most helpful way to get a sense of what is happening. But that matrix would have $V(\alpha_{is})+V(\delta_{ih})+V(\varepsilon_{ish})$ on the diagonal, with the off-diagonal elements of $Cov(\alpha_{is},\alpha_{js})+Cov(\delta_{ih},\delta_{jh})$. Sometimes both or one of these will be zero.
The standard notation you see in many textbooks is confusing because it omits the $i$ index in some of the error terms. You really have $\alpha_{is}$, not $\alpha_{s}$, where $Cov(\alpha_{is},\alpha_{js})$ can be non zero, but $Cov(\alpha_{is},\alpha_{jk})=0$. This is also why this can't be remedied by putting fixed effects/dummies for each school and field.
Here is what variance-covariance matrix for two schools and two fields (STEM and Humanities) might look like (assuming we sorted the observations by the school and then by field within each school):

Here school 1 has five teachers (2 humanities and 2 STEM), and school 2 has ten teachers (5 in each field). The regions with color are allowed to have non-zero covariances (though some elements can still be zero). The regions with zeros are assumed to be uncorrelated because those teachers are from different schools AND different fields.
