Why is this linear combination of random variables from a white noise stationary? I'm looking at the following definition of a causal AR(p) (autoregressive) model:

An AR(P) model $\phi(B)x_t= \epsilon_t$ is said to be causal if it has a stationary solution
$$x_t=\epsilon_t +\sum^{\infty}_{i=1}b_i\epsilon_{t-i}$$
where the $\sum^{\infty}_{i=1}|{b_i}|<\infty$ and $\{\epsilon_t\}$ is a white noise process.

How do we know that this solution is stationary?
 A: There are technical issues associated with the infinite sum.  Identifying and resolving them requires some attention to the definitions, so let's begin there.

*

*A white noise process $\epsilon_t$ (in discrete time indexed by the integers) when (a) its mean $E[\epsilon_t]=0$ for all $t$ and (b) its autocorrelation function $$R_\epsilon(h,t)=E[\epsilon_t\epsilon_{t+h}]$$ is finite for all integral $h$ and zero when $h\ne 0.$  See https://en.wikipedia.org/wiki/White_noise#Discrete-time_white_noise.  Accordingly, we may drop references to $t$ in the notation for the autocorrelation.


*A stationary process $X_t$ is one for which given any dimension $n\ge 1$ and any vector $(h_1,\ldots, h_n)$ of integers, the joint distribution of $(X_{t+h_1}, X_{t+h_2}, \ldots, X_{t+h_n})$ does not vary with $t.$  See https://en.wikipedia.org/wiki/Stationary_process#Definition.
(There are weaker definitions of stationarity.  You will have no trouble establishing the result for any of them by emulating the argument used here for this strong stationarity condition.)


*A countable linear combination of random variables $X_0, X_1, X_2, \ldots, X_n, \ldots$ given by coefficients $b_0, b_1, b_2, \ldots$ is the limit (if it exists) of the finite partial sums, $$\sum_{i=0}^\infty b_i X_i = \lim_{n\to\infty} \sum_{i=1}^n b_i X_i.\tag{*}$$
But limit in what sense?  The strictest sense is the pointwise limit of the random variables; but the one most germane to models is the limit in distribution.  It doesn't matter which definition we adopt, as you will soon see.

First notice that the definition of a white noise process is so general that it doesn't even guarantee a white noise process is stationary!  When we begin with a non-stationary process $\epsilon_t,$ the stated result will generally be false.  As a counterexample, let $b_i=0$ for all $i:$ it asserts $\epsilon_t + \sum b_i \epsilon_{t-i} = \epsilon_t$ is stationary, flatly contradicting the assumption.
To make any progress, then, we must assume that $\epsilon_t$ is a stationary white noise process.
With this assumption, the derivation I gave at https://stats.stackexchange.com/a/566582/919 shows that any finite linear combination of a stationary process is stationary.  Applying this result to the situation in the question with the coefficients $(1, b_1, \ldots, b_{n})$ shows

The process $X^{(n)}_t = \epsilon_t + \sum_{i=1}^n b_i \epsilon_{t-i}$ is stationary for $n=1, 2, \ldots.$

Thus, every partial sum on the right hand side of $(*)$ is a stationary process.
We need to show the limit of a sequence of stationary processes is itself stationary.  But this is trivial: the definition of stationarity means all finite $n$-variate marginal distributions are unchanged by time translation, so upon taking limits we deduce all finite $n$-variate marginal distributions of the limit are unchanged, whence (by definition of stationarity) the limiting process is stationary.  Applying this to the sequence $n\to X^{(n)}$ of processes implies
$$X_t = \lim_{n\to\infty} X_t^{(n)} = \epsilon_t + \sum_{i=1}^\infty b_i \epsilon_{t-i}$$
is stationary, QED.
