Variance of an AR and ARMA process derived from lag notation This question concerns the asymptotic variance of an $\text{ARMA}(p,q)$ process. 
Suppose that an $\text{ARMA}$ process can be rewritten as an $\text{MA}(\infty)$ process, and from this we can in principle derive the variance.  First, let's start with a stationary $\text{AR}(1)$ process, whose variance is given by:
$$\mathbb{V}(X_t) = \frac{\sigma^2}{1-\phi^2}.$$
If we rewrite the model in lag notation we have $(1-\phi L) X_t = \epsilon_t$ which we can invert to get to:
$$X_t = \frac{\epsilon_t}{1-\phi L}.$$
Now if we take the variance of both sides, this should yield:
$$\mathbb{V}(X_t) = \frac{\mathbb{V}(\epsilon_t)}{(1-\phi L )^{2}}
= \frac{\mathbb{V}(\epsilon_t)}{(1-\phi)^{2}}
\neq \frac{\mathbb{V}(\epsilon_t)}{1-\phi^2}.$$
This already confuses me, since it gives a different expression than the proper expression for the variance.  It is clear that if we re-write this as a geometric sum and take variances of the error terms (noticing that covariance of errors at different times is zero due to IID assumption) we receive the correct expression for the variance.  How is it then possible that $(1-\phi L)^2 = 1-\phi^2$ when evaluating $(1-\phi L)$ at $L=1$? 
Furthermore, in the general case of a stationary and invertible $\text{ARMA}(p,q)$ model, textbooks (like Hamilton) indicate that the variance can be derived by inverting the $\text{AR}$ polynomial, giving us a $\Psi$ polynomial (Wold decomposition) and then squaring the value of this polynomial evaluated at $L=1$ multiplied by the error variance --- i.e. $\mathbb{V}(X_t) = \sigma^2 \cdot \Psi(1)^2$.  It seems that already in the simple case of an $\text{AR}(1)$ model this argument already doesn't work, so how could it possibly work in the much more complicated case of a general $\text{ARMA}(p,q)$?
 A: Preliminary matters involving lag operators: The lag operator is a function that operates on a sequence of values to push the index $t \in \mathbb{Z}$ backward by one point in time.  Despite a notational convention that makes it look like it does, the lag operator does not operate on the actual values in the sequence, but the sequence itself.  Thus, for a sequence $\mathbf{x} \equiv \{ x_t | t \in \mathbb{Z} \}$ we get the lagged sequence:
$$L \mathbf{x} = \{ x_{t-1} | t \in \mathbb{Z} \}.$$
When we refer to $L x_t = x_{t-1}$ this is actually a shorthand notation for $L x_t \equiv (L \mathbf{x})_t = x_{t-1}$ --- i.e., the lag operator is not actually operating on the value $x_t$ even though the notation makes it look like it does.  Instead, it operates on the sequence, and then we substitute the time index afterwards.
This brings us to the issue of applying a lag operator to a constant.  Strictly speaking, the operator only operates on a sequence indexed by $t \in \mathbb{Z}$, but we generally allow it to operate on a constant by using the convention of treating the constant as a constant sequence.  That is, for any constant $\phi \in \mathbb{R}$ we suppose that there is a corresponding sequence $\boldsymbol{\phi} \equiv \{ \phi_t | t \in \mathbb{Z} \}$ with $\phi_t = \phi$ for all $t \in \mathbb{Z}$.  We then have $L \boldsymbol{\phi} = \{ \phi_{t-1} | t \in \mathbb{Z} \}$, so we get $L \phi \equiv (L \boldsymbol{\phi})_t = \phi_{t-1} = \phi$.  By treating the constant as if it were a sequence of constants, we can apply the lag operator and this moves the time index backward by one, which does not change the constant value that is the output ---i.e., we have the rule $L \phi = \phi$ for any constant $\phi \in \mathbb{R}$.

Your error using the lag polynomial: Since $L \phi = \phi$ you also have $(1-L\phi)^2 = (1-\phi)^2$, so your error does not lie in the step you have mentioned in your question.  Rather, your error lies in the implicit use of the (erroneous) step where you take:
$$\mathbb{V} \bigg( \frac{\epsilon_t}{1-\phi L} \bigg) = \frac{\mathbb{V} (\epsilon_t)}{(1-\phi L)^2}
\quad \quad \quad (\text{Erroneous equation})$$
This is part of your more general error of assuming that $\mathbb{V}(\Psi(L) \epsilon_t) = \Psi(1)^2 \mathbb{V}(\epsilon_t)$, which is not a valid equation.  By using this equation you are implicitly treating $\mathbb{V}$ as if it were a linear operator, which it is not.  The variance is a quadratic operator, and so when the lag operator shifts the error terms, this has a nonlinear effect on the variance operation.  To see what really happens, let's expand out $(1-\phi L)^{-1}$ into its geometric expansion to give:
$$\begin{aligned}
\mathbb{V} \bigg( \frac{\epsilon_t}{1-\phi L} \bigg) 
&= \mathbb{V} \bigg( \sum_{k=0}^\infty \phi^k L^k \epsilon_t \bigg) \\[6pt]
&= \mathbb{V} \bigg( \sum_{k=0}^\infty \phi^k \epsilon_{t-k} \bigg) \\[6pt]
&= \sum_{k=0}^\infty \phi^{2k} \mathbb{V}(\epsilon_{t-k}) \\[6pt]
&= \sigma^2 \sum_{k=0}^\infty \phi^{2k} \\[6pt]
&= \frac{\sigma^2}{1-\phi^2} \\[6pt]
&\neq \frac{\sigma^2}{(1-\phi L)^2}. \\[6pt]
\end{aligned}$$
As you can see, the quadratic effect of the variance operation means that the $\phi$ values end up getting squared but the lag operations do not.  This leads to a result that is different from the erroneous equation shown above, which you have implicitly used in your working.

Application to invertible ARMA models: You should now be able to see what happens when you apply the lag operator to the objects in your question.  However, even with this done correctly, several of the equations in your question do indeed appear to be incorrect.  (Since you have not shown us the sources of those equations, I have no idea whether you actually got these from textbooks or just misinterpreted other results.)  If we present an $\text{ARMA}(p,q)$ model in its $\text{MA}(\infty)$ form then we have:
$$X_t = \Psi(L) \varepsilon_t = \sum_{k=0}^\infty \psi_k \varepsilon_{t-k}.$$
Assuming that the error terms are IID with zero mean and variance $\mathbb{V}(\varepsilon_t) = \sigma^2$, we then have:
$$\begin{aligned}
\mathbb{V}(X_t) 
&= \mathbb{V} \bigg( \sum_{k=0}^\infty \psi_k \varepsilon_{t-k} \bigg) \\[6pt]
&= \sum_{k=0}^\infty \mathbb{V}( \psi_k \varepsilon_{t-k}) \\[6pt]
&= \sum_{k=0}^\infty \psi_k^2 \mathbb{V} ( \varepsilon_{t-k} ) \\[6pt]
&= \sigma^2 \sum_{k=0}^\infty \psi_k^2. \\[6pt]
\end{aligned}$$
Since $\Psi(1)^2 = (\sum_{k=0}^\infty \psi_k)^2$ it is not generally true that $\Psi(1)^2 = \sum_{k=0}^\infty \psi_k^2$.  Thus, if you have seen assertions that $\mathbb{V}(X_t) = \sigma^2 \Psi(1)^2$ then that is wrong.  However, it can be shown that $|\Psi(1)|<\infty$ is a sufficient condition for $\sum_{k=0}^\infty \psi_k^2 < \infty$ so that the variance of the process exists.
A: Hi: It's not a big deal but,  in the fourth line, you left out the lag operator $L$.
As far as your argument, I'm not a math person but I don't think squaring the denominator is legal because there's a lag operator in there. The evaluation at $L = 1$ has to be done AFTER the summation is carried out. Maybe someone on math.stackexchange could explain why that is ?
See page 78 of below for an example of where the evaluation at L = 1 is done at
the end. The link uses $B$ instead of $L$ but the meaning is the same.
https://www.stat.tamu.edu/~suhasini/teaching673/time_series.pdf
