On Fri, Oct 14, 2011 at 2:59 PM, <josef.p...@gmail.com> wrote: > On Fri, Oct 14, 2011 at 2:29 PM, <josef.p...@gmail.com> wrote: >> On Fri, Oct 14, 2011 at 2:18 PM, Alan G Isaac <alan.is...@gmail.com> wrote: >>> On 10/14/2011 1:42 PM, josef.p...@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 >>> >> > _______________________________________________ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion