def arma_estimate(X, P, Q, lag):
R = CORRELATION(X, maxlags=lag, norm='<STR_LIT>')<EOL>R0 = R[<NUM_LIT:0>]<EOL>MPQ = lag - Q + P<EOL>N = len(X)<EOL>Y = np.zeros(N-P, dtype=complex)<EOL>for K in range(<NUM_LIT:0>, MPQ):<EOL><INDENT>KPQ = K + Q - P+<NUM_LIT:1><EOL>if KPQ < <NUM_LIT:0>:<EOL><INDENT>Y[K] = R[-KPQ].conjugate()<EOL><DEDENT>if KPQ == <NUM_LIT:0>:<EOL><INDENT>Y[K] = R0<EOL><DEDENT>if KPQ > <NUM_LIT:0>:<EOL><INDENT>Y[K] = R[KPQ]<EOL><DEDENT><DEDENT>Y.resize(lag)<EOL>if P <= <NUM_LIT:4>:<EOL><INDENT>res = arcovar_marple(Y.copy(), P)    <EOL>ar_params = res[<NUM_LIT:0>]<EOL><DEDENT>else:<EOL><INDENT>res = arcovar(Y.copy(), P)    <EOL>ar_params = res[<NUM_LIT:0>]<EOL><DEDENT>Y.resize(N-P)<EOL>for k in range(P, N):<EOL><INDENT>SUM = X[k]<EOL>for j in range(<NUM_LIT:0>, P):<EOL><INDENT>SUM = SUM + ar_params[j] * X[k-j-<NUM_LIT:1>]   <EOL><DEDENT>Y[k-P] = SUM<EOL><DEDENT>ma_params, rho = ma(Y, Q, <NUM_LIT:2>*Q)     <EOL>return ar_params, ma_params, rho<EOL>
Autoregressive and moving average estimators.

    This function provides an estimate of the autoregressive
    parameters, the moving average parameters, and the driving
    white noise variance of  an ARMA(P,Q) for a complex or real data sequence.

    The parameters are estimated using three steps:

        * Estimate the AR parameters from the original data based on a least
          squares modified Yule-Walker technique,
        * Produce a residual time sequence by filtering the original data
          with a filter based on the AR parameters,
        * Estimate the MA parameters from the residual time sequence.

    :param array X: Array of data samples (length N)
    :param int P: Desired number of AR parameters
    :param int Q: Desired number of MA parameters
    :param int lag: Maximum lag to use for autocorrelation estimates

    :return:
        * A     - Array of complex P AR parameter estimates
        * B     - Array of complex Q MA parameter estimates
        * RHO   - White noise variance estimate

    .. note::
      *  lag must be >= Q (MA order)

    **dependencies**:
        * :meth:`spectrum.correlation.CORRELATION`
        * :meth:`spectrum.covar.arcovar`
        * :meth:`spectrum.arma.ma`

    .. plot::
        :width: 80%
        :include-source:

        from spectrum import arma_estimate, arma2psd, marple_data
        import pylab

        a,b, rho = arma_estimate(marple_data, 15, 15, 30)
        psd = arma2psd(A=a, B=b, rho=rho, sides='centerdc', norm=True)
        pylab.plot(10 * pylab.log10(psd))
        pylab.ylim([-50,0])

    :reference: [Marple]_