def setup_class(cls):

        nobs = 500
        ar = [1, -0.5, 0.1]
        ma = [1, 0.7]
        dist = lambda n: np.random.standard_t(3, size=n)
        np.random.seed(8659567)
        x = arma_generate_sample(ar, ma, nobs, sigma=1, distrvs=dist,
                                 burnin=500)

        mod = TArma(x)
        order = (2, 1)
        cls.res_ls = mod.fit(order=order)
        cls.res = mod.fit_mle(order=order,
                              start_params=np.r_[cls.res_ls[0], 5, 1],
                              method='nm', disp=False)

        cls.res1_table = np.array(
          [[  0.46157133,  -0.07694534,   0.70051876,  2.88693312,  0.97283396],
           [  0.04957594,   0.04345499,   0.03492473,  0.40854823,  0.05568439],
           [  9.31038915,  -1.7706905 ,  20.05795605,  7.06632146, 17.47049812],
           [  0.        ,   0.07661218,   0.        ,  0.        ,  0.        ],
           [  0.05487968,   0.04213054,   0.03102404,  0.37860956,  0.05228474],
           [  0.04649728,   0.04569133,   0.03990779,  0.44315449,  0.05996759]])

        cls.res1_conf_int = np.array([[ 0.36440426,  0.55873839],
                                   [-0.16211556,  0.00822488],
                                   [ 0.63206754,  0.76896998],
                                   [ 2.08619331,  3.68767294],
                                   [ 0.86369457,  1.08197335]])


        cls.ls_params = np.array([ 0.43393123, -0.08402678,  0.73293058])
        cls.ls_bse = np.array([ 0.0377741 ,  0.03567847,  0.02744488])

def SimulatePC1(p=50):
	m = np.zeros((nyears,p))
	for j in range(p):
		b = arma_generate_sample(ar, ma, nyears)
		m[:,j] = b - np.mean(b[500:600])
	u,s,v = svd(m,full_matrices=False)
	return(u[:,0])

def mcarma22(niter=10, nsample=1000, ar=None, ma=None, sig=0.5):
    '''run Monte Carlo for ARMA(2,2)

    DGP parameters currently hard coded
    also sample size `nsample`

    was not a self contained function, used instances from outer scope
      now corrected

    '''
    #nsample = 1000
    #ar = [1.0, 0, 0]
    if ar is None:
        ar = [1.0, -0.55, -0.1]
    #ma = [1.0, 0, 0]
    if ma is None:
        ma = [1.0,  0.3,  0.2]
    results = []
    results_bse = []
    for _ in range(niter):
        y2 = arma_generate_sample(ar,ma,nsample+1000, sig)[-nsample:]
        y2 -= y2.mean()
        arest2 = Arma(y2)
        rhohat2a, cov_x2a, infodict, mesg, ier = arest2.fit((2,2))
        results.append(rhohat2a)
        err2a = arest2.geterrors(rhohat2a)
        sige2a = np.sqrt(np.dot(err2a,err2a)/nsample)
        #print('sige2a', sige2a,
        #print('cov_x2a.shape', cov_x2a.shape
        #results_bse.append(sige2a * np.sqrt(np.diag(cov_x2a)))
        if cov_x2a is not None:
            results_bse.append(sige2a * np.sqrt(np.diag(cov_x2a)))
        else:
            results_bse.append(np.nan + np.zeros_like(rhohat2a))
    return np.r_[ar[1:], ma[1:]], np.array(results), np.array(results_bse)

def test_arma_generate_sample(dist, ar, ma):
    # Test that this generates a true ARMA process
    # (amounts to just a test that scipy.signal.lfilter does what we want)
    T = 100
    np.random.seed(1234)
    eta = dist(T)

    # rep1: from module function
    np.random.seed(1234)
    rep1 = arma_generate_sample(ar, ma, T, distrvs=dist)
    # rep2: "manually" create the ARMA process
    ar_params = -1 * np.array(ar[1:])
    ma_params = np.array(ma[1:])
    rep2 = _manual_arma_generate_sample(ar_params, ma_params, eta)
    assert_array_almost_equal(rep1, rep2, 13)

def test_arma_order_select_ic():
    # smoke test, assumes info-criteria are right
    from statsmodels.tsa.arima_process import arma_generate_sample

    arparams = np.array([.75, -.25])
    maparams = np.array([.65, .35])
    arparams = np.r_[1, -arparams]
    maparam = np.r_[1, maparams]
    nobs = 250
    np.random.seed(2014)
    y = arma_generate_sample(arparams, maparams, nobs)
    res = arma_order_select_ic(y, ic=['aic', 'bic'], trend='nc')
    # regression tests in case we change algorithm to minic in sas
    aic_x = np.array([[       np.nan,  552.7342255 ,  484.29687843],
                      [ 562.10924262,  485.5197969 ,  480.32858497],
                      [ 507.04581344,  482.91065829,  481.91926034],
                      [ 484.03995962,  482.14868032,  483.86378955],
                      [ 481.8849479 ,  483.8377379 ,  485.83756612]])
    bic_x = np.array([[       np.nan,  559.77714733,  494.86126118],
                      [ 569.15216446,  496.08417966,  494.41442864],
                      [ 517.61019619,  496.99650196,  499.52656493],
                      [ 498.12580329,  499.75598491,  504.99255506],
                      [ 499.49225249,  504.96650341,  510.48779255]])
    aic = DataFrame(aic_x, index=lrange(5), columns=lrange(3))
    bic = DataFrame(bic_x, index=lrange(5), columns=lrange(3))
    assert_almost_equal(res.aic.values, aic.values, 5)
    assert_almost_equal(res.bic.values, bic.values, 5)
    assert_equal(res.aic_min_order, (1, 2))
    assert_equal(res.bic_min_order, (1, 2))
    assert_(res.aic.index.equals(aic.index))
    assert_(res.aic.columns.equals(aic.columns))
    assert_(res.bic.index.equals(bic.index))
    assert_(res.bic.columns.equals(bic.columns))

    index = pd.date_range('2000-1-1', freq='M', periods=len(y))
    y_series = pd.Series(y, index=index)
    res_pd = arma_order_select_ic(y_series, max_ar=2, max_ma=1,
                                  ic=['aic', 'bic'], trend='nc')
    assert_almost_equal(res_pd.aic.values, aic.values[:3, :2], 5)
    assert_almost_equal(res_pd.bic.values, bic.values[:3, :2], 5)
    assert_equal(res_pd.aic_min_order, (2, 1))
    assert_equal(res_pd.bic_min_order, (1, 1))

    res = arma_order_select_ic(y, ic='aic', trend='nc')
    assert_almost_equal(res.aic.values, aic.values, 5)
    assert_(res.aic.index.equals(aic.index))
    assert_(res.aic.columns.equals(aic.columns))
    assert_equal(res.aic_min_order, (1, 2))

def test(db):

    np.random.seed(12345)
    arparams = np.array([.75, -.25])
    maparams = np.array([.65, .35])
    
    arparams = np.r_[1, -arparams]
    maparam = np.r_[1, maparams]
    nobs = 250
    y = arma_generate_sample(arparams, maparams, nobs)
    
    dates = sm.tsa.datetools.dates_from_range('1980m1', length=nobs)
    y = pandas.Series(y, index=dates)
    arma_mod = sm.tsa.ARMA(y, order=(2, 2))
    arma_res = arma_mod.fit(trend='nc', disp=-1)
    
    print(arma_res.summary())

def test_arma_generate_sample():
    # Test that this generates a true ARMA process
    # (amounts to just a test that scipy.signal.lfilter does what we want)
    T = 100;
    dists = [np.random.randn]
    for dist in dists:
        np.random.seed(1234);
        eta = dist(T);
        for ar in arlist:
            for ma in malist:
                # rep1: from module function
                np.random.seed(1234);
                rep1 = arma_generate_sample(ar, ma, T, distrvs=dist);
                # rep2: "manually" create the ARMA process
                ar_params = -1*np.array(ar[1:]);
                ma_params = np.array(ma[1:]);
                rep2 = _manual_arma_generate_sample(ar_params, ma_params, eta)
                assert_array_almost_equal(rep1, rep2, 13);

def calibrate_arma(ar=3, ma=2, data=None):
    # Fitting an ARMA model for the required data
    arma = sm.tsa.ARMA(data, (ar, ma)).fit(disp=0)
    # Capturing the ARMA params
    params = arma.params
    # Splitting the params into AR and MA params
    arparams = np.r_[1, -np.array(
        params[len(params) - (ar + ma):len(params) - ma])]
    maparams = np.r_[1, np.array(params[len(params) - ma:])]

    # Creating a callback function for generating new series
    def gen(nsample):
        # Same mean and standard deviation
        return [normalvariate(np.average(arma.resid), np.std(
            arma.resid)) for i in range(nsample)]

    # Generating new time series with the same properties of data
    return MinMaxScaler((min(data), max(data))).fit_transform(
        arma_generate_sample(arparams, maparams, len(data), distrvs=gen))

def generate_kindofgarch(nobs, ar, ma, mu=1.):
    '''simulate garch like process but not squared errors in arma
    used for initial trial but produces nice graph
    '''
    #garm1, gmam1 = [0.4], [0.2]
    #pqmax = 1
#    res = np.zeros(nobs+pqmax)
#    rvs = np.random.randn(nobs+pqmax,2)
#    for t in range(pqmax,nobs+pqmax):
#        res[i] =
    #ar = [1.0, -0.99]
    #ma = [1.0,  0.5]
    #this has the wrong distribution, should be eps**2
    #TODO: use new version tsa.arima.??? instead, has distr option
    #arest = tsa.arima.ARIMA()
    #arest = tsa.arima.ARIMA  #try class method, ARIMA needs data in constructor
    from statsmodels.tsa.arima_process import arma_generate_sample
    h = arma_generate_sample(ar,ma,nobs,0.1)
    #h = np.abs(h)
    h = (mu+h)**2
    h = np.exp(h)
    err = np.sqrt(h)*np.random.randn(nobs)
    return err, h

    _attrs = {}
    _wrap_attrs = wrap.union_dicts(tsbase.TimeSeriesResultsWrapper._wrap_attrs,
                                    _attrs)
    _methods = {}
    _wrap_methods = wrap.union_dicts(
                        tsbase.TimeSeriesResultsWrapper._wrap_methods,
                        _methods)
wrap.populate_wrapper(ARMAResultsWrapper, ARMAResults)

if __name__ == "__main__":
    import numpy as np
    import statsmodels.api as sm

    # simulate arma process
    from statsmodels.tsa.arima_process import arma_generate_sample
    y = arma_generate_sample([1., -.75],[1.,.25], nsample=1000)
    arma = ARMA(y)
    res = arma.fit(trend='nc', order=(1,1))

    np.random.seed(12345)
    y_arma22 = arma_generate_sample([1.,-.85,.35],[1,.25,-.9], nsample=1000)
    arma22 = ARMA(y_arma22)
    res22 = arma22.fit(trend = 'nc', order=(2,2))

    # test CSS
    arma22_css = ARMA(y_arma22)
    res22css = arma22_css.fit(trend='nc', order=(2,2), method='css')


    data = sm.datasets.sunspots.load()
    ar = ARMA(data.endog)

import numdifftools as ndt

import statsmodels.api as sm
from statsmodels.sandbox import tsa
from statsmodels.tsa.arma_mle import Arma  # local import
from statsmodels.tsa.arima_process import arma_generate_sample

examples = ['arma']
if 'arma' in examples:

    print "\nExample 1"
    print '----------'
    ar = [1.0, -0.8]
    ma = [1.0,  0.5]
    y1 = arma_generate_sample(ar,ma,1000,0.1)
    y1 -= y1.mean() #no mean correction/constant in estimation so far

    arma1 = Arma(y1)
    arma1.nar = 1
    arma1.nma = 1
    arma1res = arma1.fit_mle(order=(1,1), method='fmin')
    print arma1res.params

    #Warning need new instance otherwise results carry over
    arma2 = Arma(y1)
    arma2.nar = 1
    arma2.nma = 1
    res2 = arma2.fit(method='bfgs')
    print res2.params
    print res2.model.hessian(res2.params)

#
# Generate data
#
print 'Generating data...'

# generate labels
labels = np.array([1] * nq + [0] * nq + [1] * nq + [0] * nq)

# generate static features
static = np.vstack((np.random.normal(5.0, 1.0, (nq, nfeatures)),
                    np.random.normal(0.0, 1.0, (nq, nfeatures)),
                    np.random.normal(0.0, 1.0, (nq, nfeatures)),
                    np.random.normal(0.0, 1.0, (nq, nfeatures))))

# generate dynamic features
dynamic_q1 = np.vstack([ap.arma_generate_sample(arparams_pos, maparams_pos, seqlen) for i in range(nq * nseqfeatures)]).reshape(nq, nseqfeatures, seqlen) 
dynamic_q2 = np.vstack([ap.arma_generate_sample(arparams_pos, maparams_pos, seqlen) for i in range(nq * nseqfeatures)]).reshape(nq, nseqfeatures, seqlen)
dynamic_q3 = np.vstack([ap.arma_generate_sample(arparams_neg, maparams_neg, seqlen) for i in range(nq * nseqfeatures)]).reshape(nq, nseqfeatures, seqlen)
dynamic_q4 = np.vstack([ap.arma_generate_sample(arparams_pos, maparams_pos, seqlen) for i in range(nq * nseqfeatures)]).reshape(nq, nseqfeatures, seqlen)
dynamic = np.vstack((dynamic_q1, dynamic_q2, dynamic_q3, dynamic_q4))


#
# Split the dataset into training and test sets
#
print 'Splitting train and test...'

# pick samples for training and for testing
train_idx = np.random.choice(range(0, nsamples), size=np.round(nsamples * 0.7, 0), replace=False)
test_idx = list(set(range(0, nsamples)) - set(train_idx))

 x = 1+np.arange(12)
 for blocksize in [2,2,3,3,4,4,5,5,5]:
     print(blocksize, block_shuffle(x, blocksize))
 # tests of discretise
 print('\ntests of discretise')
 x = 1+np.arange(12)
 np.random.shuffle(x)
 print('x:', x)
 for k in [1,2,3,4,12]:
     print(k, discretise(x, k))
 # tests of discretise + entropy_rate
 # print('\ntests of discretise + entropy_rate')
 # x = np.random.randint(0,2**8,size=100000).astype(np.uint32)
 # xb = bytearray(x)
 # print('x[:4]: ',x[:4],'xb[:4]: ',xb[:4],'list(xb)[:4]',list(xb)[:4])
 # for method in ['lzma','bz2']:
 #     print(method,'undiscretised',entropy_rate(bytearray(x)))
 #     for discretisation_level in [2**j for j in [1,2,3,5,8,17]]:
 #         xd = discretise(x,discretisation_level)
 #         print(method,'discretised',discretisation_level,
 #         entropy_rate(bytearray(xd), method=method))
 # tests of dependogram
 print('\ntests of dependogram')
 x = arma_generate_sample([1,.5],[1],nsample=10000).astype(np.float64)
 xb = bytearray(x)
 print('x[:5]',x[:5])
 # print(dependogram(x, blocksizes=[1,2,3,100], nshuffles=100, nbins=4, plot=True))
 # dependogram(x, nshuffles=100, nbins=8, plot=True)
 # dependogram(x, nshuffles=100, nbins=6, plot=True)
 dependogram(x, nshuffles=100, nbins=2, plot=False, method='lzma')
 dependogram(x, nshuffles=100, nbins=2, plot=True, method='bz2')

    _attrs = {}
    _wrap_attrs = wrap.union_dicts(tsbase.TimeSeriesResultsWrapper._wrap_attrs, _attrs)
    _methods = {}
    _wrap_methods = wrap.union_dicts(tsbase.TimeSeriesResultsWrapper._wrap_methods, _methods)


wrap.populate_wrapper(ARMAResultsWrapper, ARMAResults)

if __name__ == "__main__":
    import numpy as np
    import statsmodels.api as sm

    # simulate arma process
    from statsmodels.tsa.arima_process import arma_generate_sample

    y = arma_generate_sample([1.0, -0.75], [1.0, 0.25], nsample=1000)
    arma = ARMA(y)
    res = arma.fit(trend="nc", order=(1, 1))

    np.random.seed(12345)
    y_arma22 = arma_generate_sample([1.0, -0.85, 0.35], [1, 0.25, -0.9], nsample=1000)
    arma22 = ARMA(y_arma22)
    res22 = arma22.fit(trend="nc", order=(2, 2))

    # test CSS
    arma22_css = ARMA(y_arma22)
    res22css = arma22_css.fit(trend="nc", order=(2, 2), method="css")

    data = sm.datasets.sunspots.load()
    ar = ARMA(data.endog)
    resar = ar.fit(trend="nc", order=(9, 0))

#Example TSA ARMA

import numpy as np
import statsmodels.api as sm


# Generate some data from an ARMA process
from statsmodels.tsa.arima_process import arma_generate_sample
arparams = np.array([.75, -.25])
maparams = np.array([.65, .35])
# The conventions of the arma_generate function require that we specify a
# 1 for the zero-lag of the AR and MA parameters and that the AR parameters
# be negated.
arparams = np.r_[1, -arparams]
maparam = np.r_[1, maparams]
nobs = 250
y = arma_generate_sample(arparams, maparams, nobs)
plt.figure()
plt.plot(y)

#Now, optionally, we can add some dates information. For this example,
# we'll use a pandas time series.
dates = sm.tsa.datetools.dates_from_range('1980m1', length=nobs)
y = Series(y, index=dates)
arma_mod = sm.tsa.ARMA(y, freq='M')
#arma_res = arma_mod.fit(order=(2,2), trend='nc', disp=-1) #fails
#old pandas 0.4.0: AttributeError: 'TimeSeries' object has no attribute 'name'
#arma_res.params

plt.show()

def generate_lstm_wins(nsamples, nfeatures, nseqfeatures, seqlen):
    
    # size of 1/4th of a dataset
    nq = nsamples / 4
    
    # dynamical feature distribution parameters (for arma)
    orders = [np.random.randint(1, 5, 4) for rep in range(nseqfeatures)]
    arparams_pos = [np.hstack((1.0, np.random.uniform(-0.1, 0.1, ord[0]))) for ord in orders]
    arparams_neg = [np.hstack((1.0, np.random.uniform(-0.1, 0.1, ord[1]))) for ord in orders]
    maparams_pos = [np.hstack((1.0, np.random.uniform(-0.5, 0.5, ord[2]))) for ord in orders]
    maparams_neg = [np.hstack((1.0, np.random.uniform(-0.5, 0.5, ord[3]))) for ord in orders]

    #
    # Generate data
    #
    print 'Generating data...'

    # generate labels
    labels = np.array([1] * nq + [0] * nq + [1] * nq + [0] * nq)

    # generate distribution parameters for the static features
    smean_pos = np.random.uniform(0, 2, nfeatures)
    ssigm_pos = [1.0] * nfeatures
    smean_neg = np.random.uniform(0, 2, nfeatures)
    ssigm_neg = [1.0] * nfeatures

    # generate static features
    static_q1 = np.vstack([np.random.normal(smean_pos[f], ssigm_pos[f], nq) for f in range(nfeatures)])
    static_q2 = np.vstack([np.random.normal(smean_neg[f], ssigm_neg[f], nq) for f in range(nfeatures)])
    static_q3 = np.vstack([np.random.normal(smean_neg[f], ssigm_neg[f], nq) for f in range(nfeatures)])
    static_q4 = np.vstack([np.random.normal(smean_neg[f], ssigm_neg[f], nq) for f in range(nfeatures)])
    static = np.vstack((static_q1.T, static_q2.T, static_q3.T, static_q4.T))

    # generate dynamic features
    dynamic_q1 = np.vstack([ap.arma_generate_sample(arparams_pos[i % nseqfeatures], maparams_pos[i % nseqfeatures], seqlen)
                           for i in range(nq * nseqfeatures)]).reshape(nq, nseqfeatures, seqlen) 
    dynamic_q2 = np.vstack([ap.arma_generate_sample(arparams_pos[i % nseqfeatures], maparams_pos[i % nseqfeatures], seqlen)
                           for i in range(nq * nseqfeatures)]).reshape(nq, nseqfeatures, seqlen)
    dynamic_q3 = np.vstack([ap.arma_generate_sample(arparams_neg[i % nseqfeatures], maparams_neg[i % nseqfeatures], seqlen)
                           for i in range(nq * nseqfeatures)]).reshape(nq, nseqfeatures, seqlen)
    dynamic_q4 = np.vstack([ap.arma_generate_sample(arparams_pos[i % nseqfeatures], maparams_pos[i % nseqfeatures], seqlen)
                           for i in range(nq * nseqfeatures)]).reshape(nq, nseqfeatures, seqlen)
    dynamic = np.vstack((dynamic_q1, dynamic_q2, dynamic_q3, dynamic_q4))


    #
    # Split the dataset into training and test sets
    #
    print 'Splitting train and test...'

    # pick samples for training and for testing
    train_idx = np.random.choice(range(0, nsamples), size=np.round(nsamples * 0.7, 0), replace=False)
    test_idx = list(set(range(0, nsamples)) - set(train_idx))

    train_static = static[train_idx, :]
    train_dynamic = dynamic[train_idx, :, :]
    test_static = static[test_idx, :]
    test_dynamic = dynamic[test_idx, :, :]
    train_labels = labels[train_idx]
    test_labels = labels[test_idx]

    return train_static, train_dynamic, test_static, test_dynamic, train_labels, test_labels

'''

from __future__ import print_function
import numpy as np
import statsmodels.api as sm
from statsmodels.tsa.arima_process import arma_generate_sample
from statsmodels.tsa.arma_mle import Arma as Arma
from statsmodels.tsa.arima_process import ARIMA as ARIMA_old
from statsmodels.sandbox.tsa.garch import Arma as Armamle_old


print("\nExample 1")
ar = [1.0,  -0.6, 0.1]
ma = [1.0,  0.5, 0.3]
nobs = 1000
y22 = arma_generate_sample(ar, ma, nobs+1000, 0.5)[-nobs:]
y22 -= y22.mean()
start_params = [0.1, 0.1, 0.1, 0.1]
start_params_lhs = [-0.1, -0.1, 0.1, 0.1]

print('truelhs', np.r_[ar[1:], ma[1:]])





###bug in current version, fixed in Skipper and 1 more
###arr[1:q,:] = params[p+k:p+k+q]  # p to p+q short params are MA coeffs
###ValueError: array dimensions are not compatible for copy
##from statsmodels.tsa.arima import ARMA as ARMA_kf
##arma22 = ARMA_kf(y22, constant=False, order=(2,2))

"""
Autoregressive Moving Average (ARMA) Model
"""
import numpy as np
import statsmodels.api as sm

# Generate some data from an ARMA process
from statsmodels.tsa.arima_process import arma_generate_sample

np.random.seed(12345)
arparams = np.array([.75, -.25])
maparams = np.array([.65, .35])

# The conventions of the arma_generate function require that we specify a
# 1 for the zero-lag of the AR and MA parameters and that the AR parameters
# be negated.
ar = np.r_[1, -arparams]
ma = np.r_[1, maparams]
nobs = 250
y = arma_generate_sample(ar, ma, nobs)

# Now, optionally, we can add some dates information. For this example,
# we'll use a pandas time series.
import pandas
dates = sm.tsa.datetools.dates_from_range('1980m1', length=nobs)
y = pandas.TimeSeries(y, index=dates)
arma_mod = sm.tsa.ARMA(y, order=(2, 2))
arma_res = arma_mod.fit(trend='nc', disp=-1)

'''

import numpy as np
from statsmodels.tsa.arima_process import arma_generate_sample
from .maketests_mlabwrap import HoldIt


if __name__ == '__main__':
    filen = 'savedrvs_tmp.py'
    np.set_printoptions(precision=14, linewidth=100)


    # check arma to return same as random.normal
    np.random.seed(10000)
    xo = arma_generate_sample([1], [1], nsample=100)
    xo2 = np.round(xo*1000).astype(int)
    np.random.seed(10000)
    rvs = np.random.normal(size=100)
    rvs2 = np.round(xo*1000).astype(int)
    assert (xo2==rvs2).all()

    nsample = 1000
    data =  HoldIt('rvsdata')

    np.random.seed(10000)
    xo = arma_generate_sample([1, -0.8, 0.5], [1], nsample=nsample)
    data.xar2 = np.round(xo*1000).astype(int)
    np.random.seed(10000)
    xo = np.random.normal(size=nsample)
    data.xnormal = np.round(xo*1000).astype(int)

try:
    import scikits.talkbox.spectral.basic as stbs
except ImportError:
    hastalkbox = False

ar = [1., -0.7]#[1,0,0,0,0,0,0,-0.7]
ma = [1., 0.3]

ar = np.convolve([1.]+[0]*50 +[-0.6], ar)
ar = np.convolve([1., -0.5]+[0]*49 +[-0.3], ar)

n_startup = 1000
nobs = 1000
# throwing away samples at beginning makes sample more "stationary"

xo = arma_generate_sample(ar,ma,n_startup+nobs)
x = xo[n_startup:]


plt.figure()
plt.plot(x)

rescale = 0

w, h = signal.freqz(ma, ar)
sd = np.abs(h)**2/np.sqrt(2*np.pi)

if np.sum(np.isnan(h)) > 0:
    # this happens with unit root or seasonal unit root'
    print('Warning: nan in frequency response h')
    h[np.isnan(h)] = 1.