[Numpy-discussion] simulate AR

Fri Oct 14 15:45:29 CDT 2011

On Fri, Oct 14, 2011 at 2:59 PM,  <josef.pktd@gmail.com> wrote:
> On Fri, Oct 14, 2011 at 2:29 PM,  <josef.pktd@gmail.com> wrote:
>> On Fri, Oct 14, 2011 at 2:18 PM, Alan G Isaac <alan.isaac@gmail.com> wrote:
>>> On 10/14/2011 1:42 PM, josef.pktd@gmail.com wrote:
>>>> If I remember correctly, signal.lfilter doesn't require stationarity,
>>>> but handling of the starting values is a bit difficult.
>>>
>>>
>>> Hmm.  Yes.
>>> AR(1) is trivial, but how do you handle higher orders?
>>
>> I would have to look for it.
>> You can invert the stationary part of the AR polynomial with the numpy
>> polynomial classes using division. The main thing is to pad with
>> enough zeros corresponding to the truncation that you want. And in the
>> old classes to watch out because the order is reversed high to low
>> instead of low to high or the other way around.
>
> I found it in the examples folder (pure numpy)
> https://github.com/statsmodels/statsmodels/blob/master/scikits/statsmodels/examples/tsa/lagpolynomial.py
>
>>>> ar = LagPolynomial([1, -0.8])
>>>> ma = LagPolynomial([1])
>>>> ma.div(ar)
> Polynomial([ 1. ,  0.8], [-1.,  1.])
>>>> ma.div(ar, maxlag=50)
> Polynomial([  1.00000000e+00,   8.00000000e-01,   6.40000000e-01,
>         5.12000000e-01,   4.09600000e-01,   3.27680000e-01,
>         2.62144000e-01,   2.09715200e-01,   1.67772160e-01,
>         1.34217728e-01,   1.07374182e-01,   8.58993459e-02,
>         6.87194767e-02,   5.49755814e-02,   4.39804651e-02,
>         3.51843721e-02,   2.81474977e-02,   2.25179981e-02,
>         1.80143985e-02,   1.44115188e-02,   1.15292150e-02,
>         9.22337204e-03,   7.37869763e-03,   5.90295810e-03,
>         4.72236648e-03,   3.77789319e-03,   3.02231455e-03,
>         2.41785164e-03,   1.93428131e-03,   1.54742505e-03,
>         1.23794004e-03,   9.90352031e-04,   7.92281625e-04,
>         6.33825300e-04,   5.07060240e-04,   4.05648192e-04,
>         3.24518554e-04,   2.59614843e-04,   2.07691874e-04,
>         1.66153499e-04,   1.32922800e-04,   1.06338240e-04,
>         8.50705917e-05,   6.80564734e-05,   5.44451787e-05,
>         4.35561430e-05,   3.48449144e-05,   2.78759315e-05,
>         2.23007452e-05], [-1.,  1.])
>
>>>> ar = LagPolynomial([1, -0.8, 0.2, 0.1, 0.1])
>>>> ma.div(ar, maxlag=50)
> Polynomial([  1.00000000e+00,   8.00000000e-01,   4.40000000e-01,
>         9.20000000e-02,  -1.94400000e-01,  -2.97920000e-01,
>        -2.52656000e-01,  -1.32300800e-01,  -6.07744000e-03,
>         7.66558080e-02,   1.01035814e-01,   7.93353139e-02,
>         3.62032515e-02,  -4.67362386e-03,  -2.90166622e-02,
>        -3.38324615e-02,  -2.44155995e-02,  -9.39695872e-03,
>         3.65046531e-03,   1.06245701e-02,   1.11508188e-02,
>         7.37039040e-03,   2.23864501e-03,  -1.86070097e-03,
>        -3.78841070e-03,  -3.61949191e-03,  -2.17570579e-03,
>        -4.51755084e-04,   8.14527351e-04,   1.32149267e-03,
>         1.15703475e-03,   6.25052041e-04,   5.50326804e-05,
>        -3.28837006e-04,  -4.52284820e-04,  -3.64068927e-04,
>        -1.73417745e-04,   1.21917720e-05,   1.26072341e-04,
>         1.52168186e-04,   1.12642678e-04,   4.58540937e-05,
>        -1.36693133e-05,  -4.65873557e-05,  -5.03856990e-05,
>        -3.42095661e-05], [-1.,  1.])

I just realized that my LagPolynomial has a filter method

>>> marep = ma.div(ar, maxlag=50)
>>> u = np.random.randn(5000)
>>> x = marep.filter(u)[1000:]
>>> import  scikits.statsmodels.api as sm
>>> sm.tsa.AR(x).fit(4, trend='nc').params
array([ 0.80183437, -0.22098967, -0.08484519, -0.12590277])
>>> #true (different convention)
>>> -np.array([1, -0.8, 0.2, 0.1, 0.1])[1:]
array([ 0.8, -0.2, -0.1, -0.1])

not bad, if the sample is large enough. I don't remember what numpy
polynomial use under the hood (maybe convolve)

>
> no guarantees, I don't remember how much I tested these things, but I
> spent a lot of time doing it 3 or 4 different ways.

Josef

>
> Josef
>
>>
>> I switched to using mostly lfilter, but I guess the statsmodels
>> sandbox (and the mailing list) still has my "playing with ARMA
>> polynomials" code.
>> (I think it might be pretty useful for teaching. I wished I had the
>> functions to calculate some examples when I learned this.)
>>
>> Josef
>>>
>>> Thanks,
>>> Alan
>>>
>>> _______________________________________________
>>> NumPy-Discussion mailing list
>>> NumPy-Discussion@scipy.org
>>> http://mail.scipy.org/mailman/listinfo/numpy-discussion
>>>
>>
>