[Numpy-discussion] polynomial fromroots
Skipper Seabold
jsseabold@gmail....
Sun Oct 10 12:38:18 CDT 2010
On Sat, Oct 9, 2010 at 10:36 PM, <josef.pktd@gmail.com> wrote:
> On Sat, Oct 9, 2010 at 10:01 PM, Charles R Harris
> <charlesr.harris@gmail.com> wrote:
>>
>>
>> On Sat, Oct 9, 2010 at 7:47 PM, <josef.pktd@gmail.com> wrote:
>>>
>>> I'm trying to see whether I can do this without reading the full manual.
>>>
>>> Is it intended that fromroots normalizes the highest order term
>>> instead of the lowest?
>>>
>>>
>>> >>> import numpy.polynomial as poly
>>>
>>> >>> p = poly.Polynomial([1, -1.88494037, 0.0178126 ])
>>> >>> p
>>> Polynomial([ 1. , -1.88494037, 0.0178126 ], [-1., 1.])
>>> >>> pr = p.roots()
>>> >>> pr
>>> array([ 0.53320748, 105.28741219])
>>> >>> poly.Polynomial.fromroots(pr)
>>> Polynomial([ 56.14003571, -105.82061967, 1. ], [-1., 1.])
>>> >>>
>>>
>>> renormalizing
>>>
>>> >>> p2 = poly.Polynomial.fromroots(pr)
>>> >>> p2/p2.coef[0]
>>> Polynomial([ 1. , -1.88494037, 0.0178126 ], [-1., 1.])
>>>
>>>
>>> this is, I think what I want to do, invert roots that are
>>> inside/outside the unit circle (whatever that means
>>>
>>> >>> pr[np.abs(pr)<1] = 1./pr[np.abs(pr)<1]
>>> >>> p3 = poly.Polynomial.fromroots(pr)
>>> >>> p3/p3.coef[0]
>>> Polynomial([ 1. , -0.54270529, 0.0050643 ], [-1., 1.])
>>>
>>
>> Wrong function ;) You defined the polynomial by its coefficients. What you
>> want to do is
>
> My coefficients are from a lag-polynomial in time series analysis
> (ARMA), and they really are the (estimated) coefficients. It is
> essentially the same as the model for scipy.signal.lfilter.
> I just need to check the roots to see whether the process is
> stationary and invertible. If one of the two lag-polynomials (moving
> average) has roots on the wrong side of the unit circle, then I can
> invert them.
>
I have just been doing
np.roots(np.r_[1,ma_params])
and the MA coefficients are invertible if
np.abs(np.roots(np.r_[1,ma_params])) < 1
That is the solution to (L**q + L**(q-1)*macoef[0] + ... +
macoef[q-1]) = 0 are inside the unit circle
Gretl, for example, uses the inverse of this so that the MA
coefficients are invertible if
(1 + macoef[0]*L + macoef[1]*L**2 + ... + macoef[q-1]*L**q) are
outside the unit circle, which uses the complex inverse
1/np.roots(np.r_[1,ma_params])
or
np.roots(np.r_[1,ma_params][::-1])
with the roots in reverse order of the original coefficients. The AR
coefficients are the same, except the parameters should be subtracted.
Ie.,
np.roots(np.r_[1,-ar_params])
should be inside the unit circle for stationarity.
Skipper
> I'm coding from memory of how this is supposed to work, so maybe I'm
> back to RTFM and RTFTB (TB=text book).
>
> (I think what I really would need is a z-transform, but I don't have
> much of an idea how to do this on a computer)
>
> Thanks, the main thing I need to do is check the convention or
> definition for the normalization. And as btw, I like that the coef are
> in increasing order
> e.g. seasonal differencing multiplied with 1 lag autoregressive
> poly.Polynomial([1.,0,0,-1])*poly.Polynomial([1,0.8])
>
> (I saw your next message: Last time I played with function
> approximation, I didn't figure out what the basis does, but it worked
> well without touching it)
>
> Josef
>
>>
>> In [1]: import numpy.polynomial as poly
>>
>> In [2]: p = poly.Polynomial.fromroots([1, -1.88494037, 0.0178126 ])
>>
>> In [3]: p
>> Out[3]: Polynomial([ 0.03357569, -1.90070346, 0.86712777, 1. ],
>> [-1., 1.])
>>
>> In [4]: p.roots()
>> Out[4]: array([-1.88494037, 0.0178126 , 1. ])
>>
>> Chuck
>>
>>
>> _______________________________________________
>> 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
>
More information about the NumPy-Discussion
mailing list