def test_acovf_fft_vs_convolution(demean, unbiased):
    np.random.seed(1)
    q = np.random.normal(size=100)

    F1 = acovf(q, demean=demean, unbiased=unbiased, fft=True)
    F2 = acovf(q, demean=demean, unbiased=unbiased, fft=False)
    assert_almost_equal(F1, F2, decimal=7)

def test_acovf2d():
    dta = sunspots.load_pandas().data
    dta.index = Index(dates_from_range('1700', '2008'))
    del dta["YEAR"]
    res = acovf(dta)
    assert_equal(res, acovf(dta.values))
    X = np.random.random((10,2))
    assert_raises(ValueError, acovf, X)

def test_acovf2d():
    dta = sunspots.load_pandas().data
    dta.index = DatetimeIndex(start='1700', end='2009', freq='A')[:309]
    del dta["YEAR"]
    res = acovf(dta)
    assert_equal(res, acovf(dta.values))
    X = np.random.random((10,2))
    assert_raises(ValueError, acovf, X)

def test_acovf_nlags_missing(acovf_data, unbiased, demean, fft, missing):
    acovf_data = acovf_data.copy()
    acovf_data[1:3] = np.nan
    full = acovf(acovf_data, unbiased=unbiased, demean=demean, fft=fft,
                 missing=missing)
    limited = acovf(acovf_data, unbiased=unbiased, demean=demean, fft=fft,
                    missing=missing, nlag=10)
    assert_allclose(full[:11], limited)

def test_acovf2d(reset_randomstate):
    dta = sunspots.load_pandas().data
    dta.index = date_range(start='1700', end='2009', freq='A')[:309]
    del dta["YEAR"]
    res = acovf(dta, fft=False)
    assert_equal(res, acovf(dta.values, fft=False))
    x = np.random.random((10, 2))
    with pytest.raises(ValueError):
        acovf(x, fft=False)

def levinson_durbin_nitime(s, order=10, isacov=False):
    '''Levinson-Durbin recursion for autoregressive processes

    '''
    #from nitime

##    if sxx is not None and type(sxx) == np.ndarray:
##        sxx_m = sxx[:order+1]
##    else:
##        sxx_m = ut.autocov(s)[:order+1]
    if isacov:
        sxx_m = s
    else:
        sxx_m = acovf(s)[:order+1]  #not tested

    phi = np.zeros((order+1, order+1), 'd')
    sig = np.zeros(order+1)
    # initial points for the recursion
    phi[1,1] = sxx_m[1]/sxx_m[0]
    sig[1] = sxx_m[0] - phi[1,1]*sxx_m[1]
    for k in xrange(2,order+1):
        phi[k,k] = (sxx_m[k]-np.dot(phi[1:k,k-1], sxx_m[1:k][::-1]))/sig[k-1]
        for j in xrange(1,k):
            phi[j,k] = phi[j,k-1] - phi[k,k]*phi[k-j,k-1]
        sig[k] = sig[k-1]*(1 - phi[k,k]**2)

    sigma_v = sig[-1]; arcoefs = phi[1:,-1]
    return sigma_v, arcoefs, pacf, phi  #return everything

def stationarity_convergence(phi):
    
        m = 20 # burn in 
        n = len(phi) - m
        if n <= 0:
                return 0
        na = max(1,int(0.1*n))
        nb = max(1,int(0.5*n))
        
        phi_a = phi[(m+n-na):(m+n)]
        phi_b= phi[(m):(m+nb)]
        phi_b= phi[(m+n-nb):(m+n)]
        
        phi_b_bar = sum(phi_b)/nb
        phi_a_bar = sum(phi_a)/na
        if len(phi_a) <= 1 or len(phi_b) <= 1:
                return 1
        
        v_a = acovf(phi_a)
        v_b = acovf(phi_b)
        
        # n gets large and na/n and nb/n stay fixed
        z_g = (phi_a_bar - phi_b_bar)/np.sqrt( v_a[0] + v_b[0] )
        return z_g

arrvs = ar_generator()
##arma = ARIMA()
##res = arma.fit(arrvs[0], 4, 0)
arma = ARIMA(arrvs[0])
res = arma.fit((4,0, 0))

print(res[0])

acf1 = acf(arrvs[0])
acovf1b = acovf(arrvs[0], unbiased=False)
acf2 = autocorr(arrvs[0])
acf2m = autocorr(arrvs[0]-arrvs[0].mean())
print(acf1[:10])
print(acovf1b[:10])
print(acf2[:10])
print(acf2m[:10])


x = arma_generate_sample([1.0, -0.8], [1.0], 500)
print(acf(x)[:20])
import statsmodels.api as sm
print(sm.regression.yule_walker(x, 10))

import matplotlib.pyplot as plt
#ax = plt.axes()

    B[:,:,1] = [[0,0],[0,0],[0,1]]
    xhat5, err5 = VARMA(x,B,C)
    #print(err5)

    #in differences
    #VARMA(np.diff(x,axis=0),B,C)


    #Note:
    # * signal correlate applies same filter to all columns if kernel.shape[1]<K
    #   e.g. signal.correlate(x0,np.ones((3,1)),'valid')
    # * if kernel.shape[1]==K, then `valid` produces a single column
    #   -> possible to run signal.correlate K times with different filters,
    #      see the following example, which replicates VAR filter
    x0 = np.column_stack([np.arange(T), 2*np.arange(T)])
    B[:,:,0] = np.ones((P,K))
    B[:,:,1] = np.ones((P,K))
    B[1,1,1] = 0
    xhat0 = VAR(x0,B)
    xcorr00 = signal.correlate(x0,B[:,:,0])#[:,0]
    xcorr01 = signal.correlate(x0,B[:,:,1])
    print(np.all(signal.correlate(x0,B[:,:,0],'valid')[:-1,0]==xhat0[P:,0]))
    print(np.all(signal.correlate(x0,B[:,:,1],'valid')[:-1,0]==xhat0[P:,1]))

    #import error
    #from movstat import acovf, acf
    from statsmodels.tsa.stattools import acovf, acf
    aav = acovf(x[:,0])
    print(aav[0] == np.var(x[:,0]))
    aac = acf(x[:,0])

def test_pandasacovf():
    s = Series(range(1, 11))
    assert_almost_equal(acovf(s), acovf(s.values))

            print "T = %i s, t = %3.2f ms"%(T,1000.*newInt)

        #try:
            bins, lc = binLightCurve(timeSpan[0],timeSpan[0]+T, times, newInt)
            lcAve = np.mean(lc)
            print "<I> = %4.2f"%lcAve
            lcVar = np.power(np.std(lc),2)
            print "sigma(I)^2 = %3.3f"%lcVar
            lcNorm = lc/lcAve
        
            #plot timestream with T total seconds binned into t=(n+1)*intTime integration time bins
            ax1.plot(bins,np.append(lcNorm,lcNorm[-1]),drawstyle='steps-post', label="%i s"%(T))

            print "Calculating auto-covariance sequence..."
            #acvs = acovf(lc,unbiased=True,demean=True) #* 1.0/(lcVar-lcAve)
            acvs = acovf(lc,unbiased=False,demean=False)

            #corr,ljb,pvalue = acf(lc,unbiased=True,qstat=True,nlags = T/newInt)
            corr,ljb,pvalue = acf(lc,unbiased=False,qstat=True,nlags = T/newInt)
            
            standalone_ljb, standalone_pvalue = acorr_ljungbox(lc)
        
            print "Min(p-value) of acf Ljung-Box test = %f"%np.min(pvalue)
            try:
                print "Min(p) of acf LB at index %i of %i"%(np.where(pvalue==np.min(pvalue))[0],len(pvalue))
                mostCorrLag = np.where(pvalue==np.min(pvalue))[0] * newInt*1000
                print "Min(p) of acf LB at lag = %4.3f ms"%mostCorrLag

            except TypeError:
                print "Min(p) of acf LB at index %i of %i"%(np.where(pvalue==np.min(pvalue))[0][0],len(pvalue))
                mostCorrLag = np.where(pvalue==np.min(pvalue))[0][0] * newInt*1000

    plt.plot(wm,sdm/sdm[0], '-', wm[maxind], sdm[maxind]/sdm[0], 'o')
else:
    plt.plot(wm, sdm, '-', wm[maxind], sdm[maxind], 'o')
plt.title('matplotlib')

if hastalkbox:
    sdp, wp = stbs.periodogram(x)
    plt.subplot(2,3,3)

    if rescale:
        plt.plot(wp,sdp/sdp[0])
    else:
        plt.plot(wp, sdp)
    plt.title('stbs.periodogram')

xacov = acovf(x, unbiased=False)
plt.subplot(2,3,4)
plt.plot(xacov)
plt.title('autocovariance')

nr = len(x)#*2/3
#xacovfft = np.fft.fft(xacov[:nr], 2*nr-1)
xacovfft = np.fft.fft(np.correlate(x,x,'full'))
#abs(xacovfft)**2 or equivalently
xacovfft = xacovfft * xacovfft.conj()

plt.subplot(2,3,5)
if rescale:
    plt.plot(xacovfft[:nr]/xacovfft[0])
else:
    plt.plot(xacovfft[:nr])

def test_acovf_warns(acovf_data):
    with pytest.warns(FutureWarning):
        acovf(acovf_data)

def test_acovf_error(acovf_data):
    with pytest.raises(ValueError):
        acovf(acovf_data, nlag=250, fft=False)

def test_acovf_nlags(acovf_data, unbiased, demean, fft, missing):
    full = acovf(acovf_data, unbiased=unbiased, demean=demean, fft=fft,
                 missing=missing)
    limited = acovf(acovf_data, unbiased=unbiased, demean=demean, fft=fft,
                    missing=missing, nlag=10)
    assert_allclose(full[:11], limited)

def test_pandasacovf():
    s = Series(lrange(1, 11))
    assert_almost_equal(acovf(s, fft=False), acovf(s.values, fft=False))