It is true that the light coming from a laser is not precisely a coherent state. One can measure the photon statistics to see that it only approximates the Poisson statistics. However, the OP is not concerned with the precise modeling of the light coming form a laser. Instead, the question is how a coherent state can come about as a result of the stimulated emission that occurs in a laser. To address this aspect, I present here a simplistic view of the process. It neglects the possibility that the excited atom remains excited and not radiate.
When the laser is switch on, it starts with the creation of a population inversion. Then one of the excited atoms decays spontaneously. The photon that is produced by the spontaneous decay stimulates other atoms to decay producing more photons. However, one can also see the subsequent stimulated decays together with the initial spontaneous decay as multiple spontaneous decays together with a normalization process.
Each spontaneous decay will effectively produce a superposition
$$ |\psi_{spon}\rangle = |\text{vac}\rangle\beta + |1\rangle\zeta , $$
where $|\text{vac}\rangle$ is the vacuum state, $|1\rangle$ is a single photon state with the correct cavity mode, and $\beta$ and $\zeta$ are coefficients. The reason for the vacuum state is not because the atom did not radiated, which I exclude here, but because some of the photons are produced in modes that would not survive in the long run. These I remove and replace by the vacuum state.
For $n$ such spontaneous decays, one gets
$$ |\psi_{n-spon}\rangle = \frac{1}{\sqrt{n!}}
\left(|\text{vac}\rangle\beta + |1\rangle\zeta\right)^n
= \frac{\beta^n}{\sqrt{n!}} \sum_{p=0}^n \frac{n!}{p! (n-p)!}
\left(|1\rangle\right)^p \frac{\zeta^p}{\beta^p} . $$
In the last expression I removed the tensor products with the vacuum states. Note that a tensor product of $p$ single photon states becomes a $p$-photon Fock state
$$\left(|1\rangle\right)^p= \sqrt{p!} |p\rangle . $$
Next, we also need to sum over $n$. For this purpose, we assume suitable $n$-dependent coefficients
$$ |\psi\rangle
= \sum_{n=0}^{\infty}\mathcal{N}_n|\psi_{n-spon}\rangle
= \mathcal{N}_0\sum_{n=0}^{\infty}\sum_{p=0}^n \frac{\beta^n}{p!(n-p)!}
\left(|1\rangle\right)^p \frac{\zeta^p}{\beta^p} , $$
where
$$ \mathcal{N}_n = \frac{\mathcal{N}_0}{\sqrt{n!}} , $$
with $\mathcal{N}_0$ being an overall normalization factor. We interchange the order of the summation and shift the one index
$$ \sum_{n=0}^{\infty}\sum_{p=0}^n F_{n,p}
= \sum_{p=0}^{\infty}\sum_{n=p}^{\infty} F_{n,p}
= \sum_{p=0}^{\infty}\sum_{q=0}^{\infty} F_{p+q,p} . $$
It leads to
$$ |\psi\rangle
= \mathcal{N}_0 \sum_{p=0}^{\infty}\sum_{q=0}^{\infty}
\frac{\zeta^p\beta^q}{p! q!} \left(|1\rangle\right)^p . $$
The sum over $q$ now combines with $\mathcal{N}_0$ to produce a new normalization constant. So, the state becomes
$$ |\psi\rangle = \mathcal{N} \sum_{p=0}^{\infty} \frac{\zeta^p}{\sqrt{p!}} |p\rangle . $$
What remains is to compute the normalization constant of the state, which would then lead to the well-known expression for the coherent state.
