A decent methodology is that employed by the automated (S)ARIMA model selector function auto.arima in "forecast" package in R. It is described in a paper by Hyndman & Khandakar "Automatic time series for forecasting: the forecast package for R" (2008). The method has been employed widely and seems to be quite popular. Its performance has been evaluated on a large set of time series and has been good, as documented in Rob J. Hyndman's blog post "R vs Autobox vs ForecastPro vs …" (see Table 1 and Table 2).
So how does it work? There are two aspects to seasonality in SARIMA modelling: seasonal differencing and seasonal AR and MA patterns. First, auto.arima uses OCSB (or optionally Canova-Hansen) test to determine whether there is a need for seasonal differencing. Second, auto.arima determines the "optimal" lag orders of SAR and SMA by doing a local search of models with different lag orders and selecting the one that minimizes the AICc (or optionally AIC or BIC).
Elaborating a little on @Richard's answer:
The model
$$
\begin{aligned}
y_t &= \gamma_0 + \gamma_2 x_t + u_t, \\
u_t &= \varphi_1 u_{t-1} + v_t.
\end{aligned}
$$
can be rearranged from noting that $u_{t} = y_t-\gamma_0 - \gamma_2 x_{t}$ so that $u_{t-1} = y_{t-1}-\gamma_0 - \gamma_2 x_{t-1}$. Plugging this into the equation for $u_t$ gives
$$
y_t-\gamma_0 - \gamma_2 x_{t}=\varphi_1(y_{t-1}-\gamma_0 - \gamma_2 x_{t-1})+v_t
$$
or
$$
y_t =\gamma_0(1-\varphi_1)+\gamma_2 x_{t}+\varphi_1y_{t-1} - \varphi_1\gamma_2 x_{t-1}+v_t
$$
Hence, the second model has a nonlinear coefficient restriction on the coefficient of the constant and, more importantly, also implicitly includes $x_{t-1}$.
Let us generate some data from the second model and estimate the restricted arima as well as unrestricted OLS model including $x_{t-1}$:
n <- 5000
x <- rnorm(n)
u <- arima.sim(list(ar=0.5),n=n)
gamma_0 <- 1
gamma_2 <- 2
y <- gamma_0 + gamma_2*x + u
The results are
> arima(y, xreg=x, order=c(1,0,0))
Call:
arima(x = y, order = c(1, 0, 0), xreg = x)
Coefficients:
ar1 intercept x
0.4880 1.0196 2.0135
s.e. 0.0123 0.0279 0.0128
and
> lm(y[2:n]~x[2:n]+y[1:(n-1)]+x[1:(n-1)])
Call:
lm(formula = y[2:n] ~ x[2:n] + y[1:(n - 1)] + x[1:(n - 1)])
Coefficients:
(Intercept) x[2:n] y[1:(n - 1)] x[1:(n - 1)]
0.5222 2.0164 0.4880 -0.9771
We see that, for example, the coefficient on $x_{t-1}$, $-0.9771$ is close to the one implied by the restricted arima model $-0.4880\cdot2.0135=-0.983$. Also, as predicted, the coefficients on $x_t$ are close in both models, and so are those on $u_{t-1}$ and $y_{t-1}$. (I suspect remaining differences to be due to different underlying algorithms in arima and lm.)
Best Answer
A seasonal ARMA model is not just a non-seasonal ARMA model of a high order, in particular of order higher than the periodicity of the data. A seasonal ARMA is defined as the product of a non-seasonal and a seasonal polynomial. Expanding the product of these polynomials we can see that the seasonal ARMA model involve additional terms that you may have overlooked at first glance.
For simplicity, let's think of an autoregressive model of order $4$ in a quarterly series.
The non-seasonal AR(4) model is defined as:
$$ \phi(L) y_t = \epsilon_t \,, \qquad \epsilon_t \sim NID(0, \sigma^2) \,, $$
where $L$ is the lag operator, thus, we can rewrite it more clearly as:
$$ (1 - \phi_1 L - \phi_2 L^2 - \phi_3 L^3 - \phi_4 L^4) y_t = \epsilon_t \,, \qquad (1) $$
or
$$ y_t = \phi_1 y_{t-1} - \phi_2 y_{t-2} - \phi_3 y_{t-3} - \phi_4 y_{t-4} y_t = \epsilon_t \,. $$
The seasonal AR(3)(1) model is defined as follows:
$$ \phi(L)\Phi(L) y_t = \epsilon_t \,, \qquad \epsilon_t \sim NID(0, \sigma^2) \,, $$
where $\phi(L)$ is the non-seasonal autoregressive polynomial $(1 - \phi_1 L - \phi_2 L^2 - \phi_3 L^3)$ and $\Phi(L)$ is the seasonal autoregressive polynomial $(1 - \phi_4 L^4)$. The product of these polynomials does not expand to the polynomial in equation (1) above. The product $\phi(L)\Phi(L)$ leads to:
$$ \phi(L)\Phi(L) = (1 - \phi_1 L - \phi_2 L^2 - \phi_3 L^3)(1 - \phi_4 L^4) y_t \\ = (1 - \phi_1 L - \phi_2 L^2 - \phi_3 L^3 - \phi_4 L^4 + \phi_1\phi_4 L^5 + \phi_2\phi_4 L^6 + \phi_3\phi_4 L^7) y_t $$
and hence it yields a polynomial different from the polynomial in equation (1).
Therefore, the state-space representation differ for these models and hence they will lead in general to different parameter estimates. You may want to try some examples with function
makeARIMAor simply compare the elementmodreturned byarimafor the models that you mention.