But is $Y = aX + cZ + e$ (as a regression model, not a math equation) also a causal model (albeit a "wrong" causal model)? If I manipulate $X$ it tells me what happens to $Y$. Doesn't it correspond to the causal graph $X \rightarrow Y, Z \rightarrow Y$ ?
It may correspond to the causal graph $X \rightarrow Y, Z \rightarrow Y$...
... But it can just as well be $Y \rightarrow X, Y \rightarrow Z$.
Statistical models are present wherever there are causal models, but those statistical models are not equivalent/identical to the causal models themselves.
A statistical model only describes the correlation, and it does not (need to) describe the causation. You can describe and fit statistical models without a description of an underlying causal model.
Or at least, certainly the statistical model alone does not tell you anything about the causation (except that there is some underlying causal mechanism, but we do not know which). In this sense it is not equivalent to a causal model.
You could see a statistical model as the shadow of a causal model.
Edit
If I write a model $Y + aX + bZ + e$ in a programming language, I could do so as a function, e.g. in Python
def model(a,X,c,Z): return a*X + b*Z + np.random.randn()```
In this case you are programming a causal model and not a statistical model. You are specifically defining the random noise as being added to the a*X+b*Z. But this does not need to be the case in order to get that there is a relation:
$$ Y|X,Z \sim N(aX + bZ, \sigma^2)$$
Consider the data below
The statistical model is that X and Y follow a bivariate distribution. But can you tell the causal model from it? Do we have $X = aY + \text{noise}$ or do we have $Y = aX +\text{noise}$ ? They can result in the same statistical distribution, but the causal models are different.
