Are all statistical models also causal models? I'm just starting to learn about causal inference methods, focused on Pearl's do-calculus.
So the point between Pearl's causal graphs and rules for manipulating causal graphs appears to be to turn a causal graph into a statistical model (e.g. a linear regression).
So you might have a causal graph such as $Z \rightarrow X, X \rightarrow M, M \rightarrow Y, Z \rightarrow Y$  (Z is a confounder of X and Y, but X also partially causes Y through M).

If my aim is to figure out the causal effect of X on Y and I just did the naive thing and setup a linear regression with $Y = aX + e$ and tried to estimate a regression coefficient, then of course I would get a biased estimate due to the presence of the confounder Z. On the other hand, if I do $Y = aX + bM + cZ + e$, then I will block the effect of X by conditioning on the mediator M. So again, knowing the causal graph will tell me to condition on Z but not on M, i.e. $Y = aX + cZ + e$ is the correct statistical model that allows me to estimate a causal effect.
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$ ? If so, then is the Pearl method just finding a transform of a causal graph into another causal graph that is easier to work with or represent as a regression?
edit: I think my causal graph analysis in this simple example is wrong, but hopefully the broader point is still clear
edit #2:
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()

So if I change the input $X$, it will cause the output of $model(...)$ to change, but I cannot do the opposite. Isn't that a causal model?


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.
