Formal definition of Simultaneity in Econometrics I'm studying Econometrics, and I'm trying to understand the endogeneity by simultaneity problem in the OLS estimation.
Given a generic OLS model:
$$y = \underline{x} \underline{\beta} + u$$ where $\underline{x}$ is the vector of explanatory variables.
My professor gave the classical, and a bit vague, explanation of simultaneity saying that it happens when an explanatory variable $x_i$ and $y$ are "codetermined", without specifying exacly what her meant by this term. I looked up the definition of simultaneity given in this wikipedia page.
In the wikipedia page it says that

Suppose that two variables are codetermined, with each affecting the other. Suppose that there are two "structural" equations,
$$\begin{align*} 
y_i &= \beta_1 x_i +\gamma_1 z_i + u_i \\ 
z_i &=  \beta_2 x_i + \gamma_2 y_i + v_i
\end{align*}
$$
[...] Assuming that $1-\gamma_1 \gamma_2 \neq 0$ and that $x_i, v_i$ are uncorrelated with $u_i$ then $$E(z_i u_i) = \dots \neq 0$$ Therefore there is endogeneity.

Now, I want to truly understand what it means for two random variables (in this case $z_i$ and $y_i$) to be codetermined. To me, this concept is highly related to the uncorrelation of the two error terms and the condition on the $\gamma$ coefficients, in fact, given any structural equation, I can always isolate an explanatory variable and write it as a linear combination of the other explanatory variables and $y_i$, but generally in this case both the condition on uncorrelation of errors and the condition on $\gamma$ will be violated, of course. So can someone explain to be, in rigorous probabilistic terms, what it means for two random variables to be codetermined (and what simultaneity is truly about)? Thanks
EDIT: I'm a math grad student, and this is the first econometric course that I'm taking, so I cannot but use a mathematical statistics approach. My question is a meta-econometric question, in the sense that the all the usual explanations I've seen (although intuitive) of this phenomenon are not precise. My question revolves around the mathematical relation between the two random variables ($y_i$ and $z_i$ in the example). For example:

*

*Does codetermination means that the random variable $z_i$ is $y_i$-measurable and vice-versa? This would imply that we can write both of them as a (deterministic) function of the other: $$z_i = f(y_i)$$ $$y_i = g(z_i)$$
 A: The simultaneity induces simultaneous equations bias. It seems you're fine with this link, and simply want a formal definition of simultaneous equations model to start with. It can be found in multiple econometrics textbooks. Here's a snippet from the chapter in Greene's 7th edition:

The system 10-42 above is the definition of simultaneous equations model. You can find all the details in the textbook, of course. You can recognize the system in your question as a special case of Eq.10-42
I suppose it's called simultaneous because you have to determine all $y_i$ at once, not one after another. However, it's more nuanced than that. For instance, seemingly unrelated regressions look quite similar to this model (except no endogenous variables on the right side), and they don't have the same endogeneity issues.
A: This is just a "static model" modelling the equilibrium
In economics, and other scientific fields, when variables are codetermined we mean that each exerts a causal effect on the other.  In such cases, there may be one or more "equilibrium" points where the mutual causal effects of the two variables hold the outcome at a fixed point.  This is common in economic processes but it also occurs in other scientific fields.  Codetermination means that each variable exerts a causal effect on the other, but this does not mean that the relationship is deterministic.  A non-deterministic relationship can occur due to exogenous causal effect operating on the variables.
What is notable about the simultaneous equations given here is that they are a "static model", meaning that they only give a description of the equilibrium, not the dynamic causal process underlying this equilibrium.  So, there is more to the process than what is described in the equations, and you are right to want to know more about the underlying causal dynamics at work.  If you wanted to model the underlying causal process leading to this equilibrium, you would use a more detailed "dynamic model" that incorporates changes over time, and the causal pressures exerted on the two variables when they are out of equilibrium.  Dynamic models of this kind are usually specified as "evolutionary" time-series models (either in discrete or continuous time) where the movement of the variables is determined by equations involving the other variable (i.e., they causally effect each other).
One of the things that some economic researchers do is to construct "dynamic models" (usually evolutionary time-series models) that describe the evolution of codetermined variables over time.  They posit a set of causal forces operating between the variables and they then investigate questions of whether or not the mutual causal effects lead to one or more equilibrium points, whether or not there is a unique equilibrium points, and whether or not the variables converge to the equilibrium(s).  They also investigate the effect of imposing external "shocks" on the variable, such as seeing how long they take to go back to equilibrium.  These are complex mathematical models, so for the most part they are used only in specialised research.  In most applied contexts, economists skip this part and use a simpler static model instead --- they directly specify simultaneous equations for a (presumed) equilibrium and they then use data to try to find the nature of the equilibrium.
A: The quantity and price of a good is a nice, simple, and hands-on example of what it means for variables to be "codetermined". I wrote this down mostly for my own benefit, but I think it is fairly relevant to post as an answer. Consider linear supply and demand
$$ Q_D = (\alpha + u) + \beta P $$
$$ Q_S = (\gamma + v) + \delta P $$
, where $u$ and $v$ are independent random variables with mean zero. In other words, the intercepts of the supply and demand curves have a constant component and a random component and will therefore shift randomly up and down. The price that we observe ($P^\star$) is determined in equilibrium $Q_D=Q_S$
$$ \alpha + \beta P^\star + u = \gamma + \delta P^\star + v $$
$$ P^\star = \frac{1}{\beta-\delta}\left(\gamma +v-\alpha-u \right) $$
If we regress $Q^\star$ on $P^\star$, the regression has two interpretations. Either we are estimating the slope of the demand curve, in which case
$$ \hat{\beta}_{OLS} \overset{p}{\to} \beta + \frac{Cov\left(P^\star, u \right)}{Var\left(P^\star\right)} = \beta + \left(\delta - \beta \right) \frac{Var(u)}{Var(u)+Var(v)} $$
or we are estimating the slope of the supply curve
$$ \hat{\delta}_{OLS} \overset{p}{\to} \delta + \frac{Cov\left(P^\star, v \right)}{Var\left(P^\star\right)} = \delta - \left(\delta - \beta \right) \frac{Var(v)}{Var(u)+Var(v)} $$
Since theory dictates that $\beta < 0$ and $\delta > 0$, we know that we have a positively biased estimator for $\beta$ and a negatively biased estimator for $\delta$. Alternatively, say that we are agnostic about what we are estimating and simply run the regression
$$ Q^\star = \lambda + \theta P^\star + \varepsilon$$
,then regardless of whether we consider $\varepsilon = u$ or $\varepsilon = v$
$$ \hat{\theta}_{OLS} \overset{p}{\to} \omega\delta + (1-\omega)\beta $$
$$ \omega = \frac{Var(u)}{Var(u)+Var(v)} $$
, the result will not correspond to any theoretical parameter of interest. Rather, a weighted average of the two parameters that we could potentially be interested in. It is evident that the bias disappears (we get something theoretically meaningful) only when either demand is constant, $Var\left(u\right)=0$, or supply is constant, $Var\left(v\right)=0$, since variation in the intercept of the one curve will trace the slope of the other.
Here is a numerical example (note that the plot below, unlike usually, has quantity on the y-axis)
alpha <- 400
beta <- -2
gamma <- -500
delta <- 10
sigmasq_u <- 1000
sigmasq_v <- 9000


If we assume $u$ and $v$ are normally distributed, and draw a random sample of 100 observations we get something like this (here, I've added a $\pm1$ standard deviation band to the curves. For teaching purposes, we could also plot just a few observations and the equilibria that generated them).

In this example, the bias we should get from simultaneity is
bias_beta <- (delta-beta)*(sigmasq_u/(sigmasq_u + sigmasq_v))
bias_delta <- -(delta-beta)*(sigmasq_v/(sigmasq_u + sigmasq_v))

> bias_beta
[1] 1.2
> bias_delta
[1] -10.8
> beta + bias_beta
[1] -0.8
> delta + bias_delta
[1] -0.8

We can simulate to confirm
n <- 100
nsim <- 10000
res <- list()

for(i in 1:nsim) {

  #print(i)
  v <- rnorm(n, 0, sqrt(sigmasq_v))
  u <- rnorm(n, 0, sqrt(sigmasq_u))
  p <- (1/(beta-delta))*(gamma + v - alpha - u)
  q <- alpha + beta*p + u # q <- gamma + delta*p + v

  res[[i]] <- coef(lm(q ~ p))

}

res <- as.data.frame(do.call("rbind", res))

> mean(res$p)
[1] -0.8003049

