ARMA model vs. Linear Filter
Based on your exploration of the topics I can make the following explanation, mostly following the approach in Statistical Digital Signal Processing and Modeling, HAYES.

Linear Time Invariant (LTI) filtering is associated with the following LCCDE (Linear Constant Coefficient Difference Equation) with zero initial conditions and assuming a causal solution; i.e., interested in the solution for $n \geq 0$ in general.

$$ \sum_{k=0}^{N} a_k y[n-k] = \sum_{k=0}^{M} b_k x[n-k]  \tag{1} $$

where $a_k$ and $b_k$ are the coefficients of the LCCDE (mostly $a_0 = 1$), and $x[n]$ is the input, while $y[n]$ is the output of the filter (aka solution of LCCDE). 

The filter is associated with an impulse response $h[n]$ (coefficients of the LTI filter) which is the solution of Eq(1) when $x[n] = \delta[n]$; $x[0] = 1$ and $x[n] = 0$ for all $n&gt; 0$. 

$$ x[n] \longrightarrow \boxed{ h[n] } \longrightarrow y[n] $$

Impulse response can be used to compute the solution $y[n]$ of the LCCDE in Eq(1) for a given input $x[n]$, through what's called a convolution sum :

$$ y[n] = h[n] \star x[n] = \sum_{k=-\infty}^{\infty} h[k] x[n-k] \tag{2} $$

Z-transform of the impulse response is the following and is useful as an algebraic transform representation of the LTI filter system.

$$ H(z) = \frac{ \sum_{k=0}^{M} b_k z^{-k} }{\sum_{k=0}^{M} a_k z^{-k}} \tag{3}$$

From which the output $y[n]$ can also be expressed as
$$ Y(z) = H(z) X(z) \tag{4} $$

which is nothing but Z-transform applied to Eq(2).

ARMA-(p,q) modeling :

Consider an LTI filter with LCCDE coefficients $a_k$ and $b_k$, as given in Eq(1), then if the input of this filter is a white-noise $v[n]$ with unit variance, then the output of this filter is called (by definition) Auto-Regressive Moving Average process with orders (p,q), abb as ARMA-(p,q), the associated LCCDE is written like 

$$ \sum_{k=0}^{p} a_k x[n-k] = \sum_{k=0}^{q} b_k v[n-k]  \tag{5} $$

Note that with the filtering convention in Eq(1), $x[n]$ was the input and $y[n]$ was the output, but withing the ARMA-(p,q) convention in Eq(5), $v[n]$ is the input and $x[n]$ (the ARMA process) is the output. The AR order is $p$ (associated with the lefthand side, recursive part, of the LCCDE) and the MA order is $q$ (associated with the right hand side of the LCCDE).

Now given some WSS RP $x[n]$, then you can model this process as an ARMA-(p,q) by finding the set of coefficents $a_k$ and $b_k$ according to Eq(5), which would minimize the mean square error between the modeling ARMA-(p,q) filter output $\hat{x}[n]$ and the given process $x[n]$ to be modeled. The solution for $a_k$ are generated by Modified Yule-Walker Equations and $b_k$ are found in a more complex manner following MA modeling techniques etc.

Linear Prediction :

Given a WSS random process $x[n]$, we want to estimate its current value $x[n]$ by  using its $p$ past values  $x[n-1],x[n-2],...,x[n-p]$ via a linear combination :

$$ \hat{x}[n] = \sum_{k=1}^{p} c_k x[n-k] $$

I have used these $c_k$ as the linear prediction coefficients; it can be shown that these coefficients are related to all-pole signal model Prony coefficients, or an ARMA-(p,0) model coefficients as $c_k = -a_k$.

Rewriting Eq(5) with $q=0$,and $a_0 = 1$, as 

$$ x[n] = - \sum_{k=1}^p a_k x[n-k] + b_0 v[n] $$

with a suitable $b_0$, reveals that, if the predicted process $x[n]$ was generated via an ARMA-(p,0) model using Eq(5), then the linear prediction coefficients $c_k$ which  minimize the prediction error will be $c_k= -a_k$ of Eq(5).