Lag operator; particular solution to ARMA as a MA$(\infty)$ process Let's take an AR(1) model. I am comfortable with the fact that $Ly_t$ means that the lag operator operates on the process ${y_t}$ by lagging it by one period. What I am a bit less comfortable with  to write a lag operator by itself without a process to operate on. In the denominator we have $1-a_1L$. It isn't immediately clear which process the lag operator is operating on, but upon a moment's reflection, I guess we realize we can multiply both sides by the denominator and I guess we can figure out that the lag is operating on $y_t$?
$$y_t=a_0+a_1Ly_t+\epsilon_t$$
$$y_t=\frac{a_0+\epsilon_t}{1-a_1L}$$
Now let's switch to the ARMA(p,q) model.
$$y_t=a_0+\sum_{i=1}^{p}a_iy_{t-i}+\sum_{i=0}^{q}\beta_i\epsilon_{t-i}$$
$$y_t-\sum_{i=1}^{p}a_iy_{t-i}=a_0+\sum_{i=0}^{q}\beta_i\epsilon_{t-i}$$
$$(1-\sum_{i=1}^{p}a_iL^{i})y_t=a_0+\sum_{i=0}^{q}\beta_i\epsilon_{t-i}$$
Because I am struggling to wrap my mind around the lag operator I cannot quite convince myself of what Enders says (page 51) "The important point to recognize is that the expansion (of the equation below) will yield an $MA(\infty)$ process". I am looking for someone to help me see that the expansion of $y_t$ expression below will yield an $MA(\infty)$ process and help me understand what it means for a lag operator to be written by itself, without something to operate on. What is in the denominator of $y_t$? How do I think about it?
$$y_t=\frac{a_0+\sum_{i=0}^{q}\beta_i\epsilon_{t-i}}{1-\sum_{i=1}^{p}a_iL^{i}}$$
 A: A brief overview of the lag operator in time-series analysis
In the answer below I will give a brief overview of the lag operator, showing its definition and properties, and how it applies in the context of time-series analysis.  For simplicity, I will show you how it works when you are dealing with time-series that are unbounded sequences of real random variables, though it is possible to alter or generalise this analysis when dealing with time-series on other spaces (e.g., complex time-series).  In any of these cases the essential working is the same, so I will take real time-series as a suitable expository context.

Defining the lag operator: The lag operator operates on a sequence of values to shift the time index by one unit.  To facilitate analysis of this operator, suppose we let $\mathbf{x} = (x_t|t \in \mathbb{Z})$ denote a sequence of real numbers and let $\mathscr{S} \equiv \mathbb{R}^\mathbb{Z}$ denote the set of all possible real sequences of this kind.$^\dagger$  For clarity, we will also define the operators $[ \ \ \ ]_i$ (for all indices $i \in \mathbb{Z}$) to recover the $i$th element from a sequence (e.g., we have $[\mathbf{x}]_i = x_i$).
With this basic apparatus set up, the lag operator is a function $L: \mathscr{S} \rightarrow \mathscr{S}$ defined by:
$$L(\mathbf{x}) = (x_{t-1}|t \in \mathbb{Z})
\quad \quad \quad \text{for all } \mathbf{x} \in \mathscr{S}.
\quad \quad \quad$$
Equivalently, the lag operator can be defined by the scalar requirement that:
$$\quad \quad \quad \ [L(\mathbf{x})]_t = [\mathbf{x}]_{t-1}
\quad \quad \quad \text{for all } \mathbf{x} \in \mathscr{S} \text{ and } t \in \mathbb{Z}.$$
By convention, we use the succinct notational equivalence $L \mathbf{x} \equiv L(\mathbf{x})$, which allows us to show the operator operating on time-series without having to put the argument series in brackets each time (though we may still need to use brackets in cases when there is a need to avoid ambiguity of the argument).  Consequently, if $\mathbf{x}$ is a sequence of real numbers then $L \mathbf{x}$.  (It is also worth noting that the lag operator sometimes goes by the name "backshift operator" and is denoted in this case by $B$.)
As you can see, the lag operator is well-defined as a function operating in the "function space" $\mathscr{S}$.  Since it is a well-defined function, it is legitimate to refer to this object without attaching an argument value on which it operates.  If you refer to it without an argument then it is a function operating in function-space.  If you refer to it with an argument attached then you are referring to the output of the function for that input argument.

A common misconception: For practitioners having difficulty understanding the lag operator in the context of time-series analysis, a common misconception is to read the lag operator as if it is mapping scalars to scalars and then suppose that the time-series is built up of a sequence of such scalars (i.e., scalars that have been operated on by the lag operator).  For example, such a practitioner might incorrectly believe that:
$$L([\mathbf{x}]_t) = [\mathbf{x}]_{t-1}.
\quad \quad \quad \text{(Erroneous equation)}$$
This misconception leads some practitioners to wonder how the lag operator can operate on identical outcomes of values in a time-series to map them to different lagged values.  That is, they note that we can have a situation where:
$$[\mathbf{x}]_t = [\mathbf{x}]_{k} 
\quad \quad \quad \quad \quad 
[\mathbf{x}]_{t-1} \neq [\mathbf{x}]_{k-1},$$
for some time-indices $t \neq k$, which would then seem to imply the contradiction:
$$L([\mathbf{x}]_t) = [\mathbf{x}]_{t-1} \neq [\mathbf{x}]_{k-1} = L([\mathbf{x}]_{k}) = L([\mathbf{x}]_t)
\quad \quad \quad \text{(Erroneous equation)}$$
The reason this does not occur is that the lag operator is mapping sequences to sequences not scalars to scalars.  So long as you bear in mind that the lag operator operates on the entire sequence, you will be aware that it can "remember" the time-placement of things and so it can map a sequence in such a way that the same scalar values in a sequence are "mapped" to different scalars in the output sequence.

Algebraic transformation of the lag operator: In order to understand the other aspect of your question, suppose we have a constant $a \in \mathbb{R}$ and we refer to an object like $1-aL$.  By convention, what we mean here is that this is a new function that is formed by taking an affine transformation of the lag operator according to the stated formula operating on the sequence (with the "one" actually referring to the identity operator in function-space).  Formally, this object is a function $(1-aL): \mathscr{S} \rightarrow \mathscr{S}$ defined by:
$$(1-aL)(\mathbf{x}) = (x_t - a x_{t-1}|t \in \mathbb{Z})
\quad \quad \quad \text{for all } \mathbf{x} \in \mathscr{S}.
\quad \quad \quad$$
Observe that this function satisfies the equivalence:
$$[(1-aL)(\mathbf{x})]_t = [\mathbf{x}]_t - a [L(\mathbf{x})_t] = [\mathbf{x}]_t - a [\mathbf{x}]_{t-1},$$
which means that the output is invariant to whether the transform is applied to the lag operator (according to our convention of what this means), or applied to the output of the lag operator.  As before, by convention we often use the succinct notational equivalence $(1-aL)\mathbf{x} \equiv (1-aL)(\mathbf{x})$ to avoid excess bracketing.
This basic idea can be extended to other transformations that operate on the entire sequence of numbers, or any scalar function that is "extended" to act as if it were a transformation operating on a sequence.  In particular, if we let $T: \mathscr{S} \rightarrow \mathscr{S}$ be an arbitrary transformation operating on the sequence then we can define the compositional operator $(TL): \mathscr{S} \rightarrow \mathscr{S}$ by:
$$(TL)(\mathbf{x}) = T(L(\mathbf{x}))
\quad \quad \quad \text{for all } \mathbf{x} \in \mathscr{S}.
\quad \quad \quad$$
In some cases (such as the affine transformation above) we appear to be using a scalar function $T_*: \mathbb{R} \rightarrow \mathbb{R}$, but we are really referring to its sequence-extension, defined by $T(\mathbf{x}) = (T_*(x_t) | t \in \mathbb{Z})$.  In this case we then have $[(TL)(\mathbf{x})]_t 
= [T(L(\mathbf{x}))]_t 
= [(T_*(x_{t-1}) | t \in \mathbb{Z})]_t 
= T_*(x_{t-1})$ for all $t \in \mathbb{Z}$, which means that the compositional operator is essentially just applying the transformation on each scalar value after the lagging operation.  (This is what is happening in the affine transformation example used above.)

Applying the lag operator to a random time-series: Our interest in the lag operator is in its use in the context of statistical time-series analysis, where we are dealing with random sequences.  In this context, our random time-series is a random sequence $\mathbf{X} : \Omega \rightarrow \mathscr{S}$ that maps each outcome in the sample space to a real sequence of values.  The lag operator can be applied to this object to yield a new random time-series $(L \mathbf{X}) : \Omega \rightarrow \mathscr{S}$.  This latter random sequence is defined by the composition:
$$(L \mathbf{X})(\omega) = L(\mathbf{X}(\omega))
\quad \quad \quad \text{for all } \omega \in \Omega.$$
In particular, it satisfies the scalar requirement that:
$$\quad \quad \quad \ [L(\mathbf{X}(\omega))]_t = [\mathbf{X}(\omega)]_{t-1}
\quad \quad \quad \text{for all } t \in \mathbb{Z} \text{ and } \omega \in \Omega.$$
From this result we can see that the lagged (random) time-series is essentially the same sequence of random variables as the original time-series, only the time-index is displaced by one in each case.  Observe in this analysis that the lag operator maps sequences to sequences and the (random) lagged time-series is defined by a composition of an initial random sequence and this sequence-to-sequence mapping.
As before, it is important to understand that the lag operator is a mapping of sequences to sequences and does not operate on individual scalars (i.e., individual random variables in the time-series).  This means that it is possible for the lag operator to "remember" the time-placement of random variables in the time-series.  It can therefore map the sequence in such a way that identical random variables in a time-series are "mapped" to different output random variables.
"Inverting" polynomials involving the lag operator: In time-series analysis we often form compositional operators that are polynomial functions of the lag operator (called "characteristic polynomials").  This is particular useful in framing ARMA models, since it allows these models to be stated and manipulated in a succinct form.  When we "invert" these characteristic polynomials, what we are really doing is finding other sequence-to-sequence operators that when multiplied with the characteristic polynomials yield the identity operator.  Fortunately, the mathematical rules for doing this turn out to be very similar to inversion of scalar polynomials, with the lag operator being treated as if it were the number one.
In this related answer I give a hueristic demonstration that the infinite moving-average form of an inverted auto-regression characteristic polynomial is an application arises from the rules for an infinite geometric sum.  If you want to look at the details when operating in function-space, you can find formal analysis for the lag operator in Kasparis 2016.

$^\dagger$ We can alter our analysis to deal with the lag operator operating on time-series that are not real numbers (e.g., complex numbers or some other kind of entity) but here we will examine real time-series.
A: Hi ColorStatistics: This isn't a mathematical proof but you can use the lag operator as if it was a number. (this is proven in functional analysis and there's a proof if ot somewhere on the net that I'll look for after I write this).
So, suppose you had $\frac{1}{1 - \rho}$ and $\rho$ was less than 1.0.
Then, using the formula for infinite geometric series, the expression
can be re-written as
$1 + \rho + \rho^2 + \rho^3 + \ldots   \rho^{\infty}$ and converges since
$\rho \lt 1.0$.
The proof of above is more straightforward than the proof that one can use
L as a number and should be in any calculus or pre-calc text.
Now, suppose you instead have an expression $\frac{\epsilon_t}{(1 - \rho L)}$ and $\rho L$ was less than 1.0 and L denotes the lag operator.
Then, you can think of $\rho L$ as the number that is acting on $\epsilon_t$.
So, since $\frac{1}{(1 - \rho L)} = $
$1 + \rho L + (\rho L)^2 + (\rho L) ^3 + \ldots   (\rho L)^{\infty}$,
one can just multiply this by $\epsilon_t$ which results in
$\epsilon_t + \rho \epsilon_{t-1} +  \rho^2 \epsilon_{t-2} + \ldots +
\rho^{\infty} \epsilon_{t-\infty}$
So, the expression above is the $MA(\infty)$ which is equivalent to the AR(1).
Another more intuitive way to think about is the following:
The AR(1) is : $y_{t} = \rho y_{t-1} + \epsilon_t$.
So, the model is such that the current response is always a combination of
some fraction, $\rho$ of the previous response + a new error term.
Note that, if one starts from the very, very beginning of the process and $y_{0} = 0$, then really what you are doing is taking an exponentially smoothed value of the error term where the exponentially smoothed weights get less and less for values of $\epsilon_{t-i}$ that are further in the past. So, the $MA(\infty)$ representation is saying: "Forget about the previous response. What one is  really doing in the AR(1) is taking new responses ( the $\epsilon_{t-i})$ and exponentially smoothing them over time
with decreasing weights.
Note that in the case of the $AR(p)$, the argument is more complex but the concept is the same. If I find the paper, I'll let you know in a comment.
