Fixed effects vs the dummy variables themselves: structural vs practical equation I have a question about if there is a substantive difference between a fixed effect and the way we estimate them (e.g., dummy variables). Are the estimated dummy variables the fixed effect, or do they simply absorb the fixed effect (and other variables invariant across the other dimensions of the data)?
I have seen fixed effect models written various ways, and I am wondering if this reflects the structural versus estimated model. For instance, one can write their fixed effects model as:
$y_{i,t} = \beta x_{i,t} + \mu_i + error$
where $\mu_i$ is the individual level fixed effect. So here, is this a 'structural model' where I can stack all $i$ observations, and have $\mu = [\mu_1,\mu_1,...\mu_n, \mu_n]'$ be one single vector which represents the fixed static $Y$ associated for each $i$?
But when we estimate this using say least squares dummy variables (LSDV), is it more accurate to write this equation as:
$y_{i,t} = \beta x_{i,t} + \sum_{i=1}^n d_i\lambda_i + error$,
where $d_1  = 1$ if $i = 1$ for all observations?
So in this specification, we aren't actually estimating the true fixed effect, we are instead using dummy variables to estimate the $\hat{\lambda_i}$'s as a way to absorb time invariant components of $y$ for each $i$? And if this understanding is correct, then is the second equation the more accurate way to write down what we are estimating in practice?
 A: 
Are the estimated dummy variables the fixed effect, or do they simply absorb the fixed effect (and other variables invariant across the other dimensions of the data)?

To be clear, estimating your equation via least squares dummy variables (LSDV) is algebraically equivalent to estimation in deviations from means. Put differently, including indicator variables for all $N-1$ entities in your panel produces mathematically equivalent estimates of $\beta$ to those where you run ordinary least squares on the 'time demeaned' data.
Furthermore, the fixed effects do not absorb variables invariant across all dimensions. In your equation, $\mu_{i}$ is $i$-subscripted and denotes cross-sections of individuals; these individuals are also observed over time. But note that $\mu_{i}$ is devoid of a $t$-subscript, and so all of the $i$-specific variation is across individuals. By restricting attention to variation within-individuals, we adjust for all time invariant confounders, both observed and unobserved.
Suppose you wanted to investigate the effect of unemployment on vehicle thefts at the city level. You sample five cities and observe general rates of auto larceny over two time periods (i.e., 2019 and 2020). Your equation would look something like the following, where I specified each city effect explicitly:
$$
\text{Larceny}_{it} = \alpha + \gamma \text{Unemployment}_{it} + \mu_{2} \text{City}_2 + \mu_{3}\text{City}_3 + \mu_{4}\text{City}_4 + \mu_{5}\text{City}_5 + \tau_{t}d2020_{t} + u_{it},
$$
where $\text{Unemployment}_{it}$ is the unemployment rate in city $i$ in time period $t$. This model includes four individual $\text{City}_{i}$ variables (i.e., dummies for all cities). Note that $\text{City}_{1}$ is dropped to avoid collinearity; it is absorbed into the intercept. The variable $d2020_{t}$ is a time dummy for the second time period.
To see why $\mu_{i}$ is the city effect, suppose you winnowed down your data frame to examine observations from the fifth city, which is unique with respect to many of its stable attributes (e.g., geography, land area, number of expressways, etc.). The equation for your city of interest simplifies to the following:
\begin{align*}
\text{Larceny}_{5t} &= \alpha + \gamma \text{Unemployment}_{5t} + \mu_{2}(0) + \mu_{3}(0) + \mu_{4}(0) + \mu_{5}(1) + \tau_{t}d2020_{t} + u_{5t} \\
\text{Larceny}_{5t} &= \alpha + \gamma \text{Unemployment}_{5t} + \mu_{5} + \tau_{t}d2020_{t} + u_{5t}
\end{align*}
where $\mu_{5}$ is the effect for the fifth city. I often like to refer to the coefficients associated with the $i$ dummies as city-specific intercepts. In simple terms, each city $i$ has it own unique starting point. You may observe variation in auto larceny rates—but it will stay constant over time for any city. In other words, the estimate of the slope coefficient (i.e., $\gamma$) is the same across all cities, but each city has its own fixed origin.
In this fictitious example, $\mu_{5}$ represents the time-constant attributes of the fifth city. For example, the geographic attributes of the fifth city represent more or less "fixed" aspects of that city. If you differenced or demeaned this equation, your city effect will be removed. Again, to see why this works, for each $i$ suppose we calculated the average of a jurisdiction's square mileage over time $t$; the average of this time-invariant variable is that time-invariant variable, so you're differencing out any fixed characteristics of that city. In fact, any variable that does not exhibit any variation over time can be safely dropped from your analysis. Your model is using the time variation within a city to identify the effect of unemployment on auto theft.
The equation is more accurately expressed in the following manner without the 'global' intercept:
$$
\text{Larceny}_{it} = \mu_{i} + \tau_t + \gamma \text{Unemployment}_{it} + u_{it},
$$
where $\mu_{i}$ denotes city fixed effects. Try 'filtering out' all observations except those from the first city and it should become even more obvious why I dropped $\alpha$. Again, $\mu_{i}$ represents a "fixed" intercept for each $i$. It's indicative of all the stable characteristics particular to each city. In practice, we often worry this unobserved city-specific heterogeneity might influence the primary variable of interest (i.e., $\text{Unemployment}_{it}$). Estimating this equation via fixed effects will remove all time-invariant confounding.
I also included year fixed effects, represented by $\tau_{t}$. It is modeling the common shocks to all $i$ in each year $t$. It could generalize to any number of years.

And if this understanding is correct, then is the second equation the more accurate way to write down what we are estimating in practice?

Your second specification uses summation notation, which denotes the estimation of dummy variables for each cross-sectional unit. I wouldn't say it's a more "accurate" estimation method though. Rather, I would argue it is more notationally explicit, as a variable (i.e., $d_{i}$) is appended to your parameter (i.e., $\lambda_{i}$). But remember, you wouldn't want to use the LSDV estimator in settings with large $N$. Suppose in your example you acquired data on 50,000 individuals observed over 20 years. Do you really want to compromise your computer storage by trying to estimate 49,999 fixed effects? Demeaning is a useful trick to overcome this.
Your post also alluded to the notion of the "true" fixed effect. But we don't know what the true fixed effect is. To be clear, $\mu_{i}$ is, in my example, an unknown parameter (intercept) for each $i$; it varies across $i$ but not over $t$. This approach adjusts for all time-constant factors, even those you haven't even considered! Thus, do not go out of your way trying to measure all entity-specific attributes invariant to time. Assuming the unobserved factors do not change over time, then the variation observed in your dependent variable is likely due to something other than these "fixed" characteristics.
A: Fixed effects
An example from this question  for a 2x2 factor model (treatment field or hedge, and site HMB or RF)
$$y_i = \beta_1 + \beta_2 \text{ treatment$_i$} + \beta_3 \text{ site$_i$} + \epsilon_i$$
effectively becomes
$$y_i = 
\left\{
  \begin{array}{@{}ll@{}}
    \beta_1, & \text{if 'treatment$_i$ = field' and 'site$_i$ = HMB'}\\
    \beta_1+\beta_2, & \text{if 'treatment$_i$ = Hedge' and 'site$_i$ = HMB'}\\
    \beta_1+\beta_3, & \text{if 'treatment$_i$ = field' and 'site$_i$ = RF'}\\
    \beta_1+\beta_2+\beta_3, & \text{if 'treatment$_i$ = Hedge' and 'site$_i$ = RF'}\\
  \end{array}\right.
$$
So indeed the vector equation with multiplication
$$\begin{bmatrix} y_1 \\ y_2 \\ \vdots  \\ y_{n-1}\\ y_{n}\end{bmatrix} = \beta_1 \cdot \begin{bmatrix} 1 \\ 1 \\ \vdots  \\ 1\\ 1\end{bmatrix} + \beta_2 \cdot \begin{bmatrix} 1 \\ 1 \\ \vdots  \\ 0\\ 0\end{bmatrix} + \beta_3 \cdot \begin{bmatrix} 0 \\ 0 \\ \vdots  \\ 1\\ 1\end{bmatrix} + \begin{bmatrix}  \epsilon_1 \\ \epsilon_2 \\ \vdots  \\ \epsilon_{n-1}\\ \epsilon_{n}\end{bmatrix}$$
is effectively also a model with the $\beta_i$ absorbed into the vectors
$$\begin{bmatrix} y_1 \\ y_2 \\ \vdots  \\ y_{n-1}\\ y_{n}\end{bmatrix} =   \begin{bmatrix} \beta_1 \\ \beta_1 \\ \vdots  \\ \beta_1\\ \beta_1\end{bmatrix} + \begin{bmatrix} \beta_2 \\ \beta_2 \\ \vdots  \\ 0\\ 0\end{bmatrix} + \begin{bmatrix} 0 \\ 0 \\ \vdots  \\ \beta_3 \\ \beta_3\end{bmatrix} + \begin{bmatrix}  \epsilon_1 \\ \epsilon_2 \\ \vdots  \\ \epsilon_{n-1}\\ \epsilon_{n}\end{bmatrix}$$
The former with multiplication is more useful. It allows you to write the matrix equation $y =  X \cdot \beta + \epsilon$
$$\begin{bmatrix} y_1 \\ y_2 \\ \vdots  \\ y_{n-1}\\ y_{n}\end{bmatrix} =   \begin{bmatrix} 1 & 1 & 0 \\ 1 & 1 & 0 \\ \vdots & \vdots & \vdots \\ 1&0&1\\ 1&0&1\end{bmatrix} \cdot \begin{bmatrix}\beta_1 \\ \beta_2 \\ \beta_3 \end{bmatrix}  + \begin{bmatrix}  \epsilon_1 \\ \epsilon_2 \\ \vdots  \\ \epsilon_{n-1}\\ \epsilon_{n}\end{bmatrix}$$
and solve the (least squares) estimate as the product of matrices (which are made up out of those dummies variables) $\beta = (X^TX)^{-1}X^T y$
Mixed effects
If the $\mu_i$ are random effects then you do not (need to) compute them (that is, you can compute $\beta$ without explicitly computing the $\mu_i$). The $\mu_i$ can be considered as a part of the error term.
I guess that it is mostly for random effects that you see the notation $y_{i,t} = \beta x_{i,t} + \mu_i + error_{i,t}$
