[Numpy-discussion] simulate AR

[josef.pktd at gmail.com]

On Fri, Oct 14, 2011 at 2:29 PM,  <josef.pktd at gmail.com> wrote:
> On Fri, Oct 14, 2011 at 2:18 PM, Alan G Isaac <alan.isaac at gmail.com> wrote:
>> On 10/14/2011 1:42 PM, josef.pktd at 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.])

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

>
> 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 at scipy.org
>> http://mail.scipy.org/mailman/listinfo/numpy-discussion
>>
>