[Numpy-discussion] polynomial fromroots
josef.pktd@gmai...
josef.pktd@gmai...
Sun Oct 10 13:19:02 CDT 2010
On Sun, Oct 10, 2010 at 1:38 PM, Skipper Seabold <jsseabold@gmail.com> wrote:
> 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.
Good, this confirms the differences in convention z, or 1/z (and why I
never remember if the roots are supposed to be inside or outside the
unit circle)
>>> arpoly.coef
array([ 1. , -0.66073887, 0.30417894])
>>> arpoly.roots()
array([ 1.08610225-1.45186791j, 1.08610225+1.45186791j])
>>> np.roots(arpoly.coef)
array([ 0.33036944+0.44162765j, 0.33036944-0.44162765j])
e.g. dividing polynomials
>>> poly.Polynomial([1,0.5,0,0,0, 0 , 0,0,0,0,0][::-1], [-1., 1.])/poly.Polynomial([1, -0.8][::-1])
Polynomial([ 0.21810381, 0.27262976, 0.3407872 , 0.425984 , 0.53248 ,
0.6656 , 0.832 , 1.04 , 1.3 , 1.
], [-1., 1.])
>>> import scikits.statsmodels.sandbox.tsa.arima as tsat
>>> tsat.arma2ma([1,-0.8],[1, 0.5],10)
array([ 1. , 1.3 , 1.04 , 0.832 , 0.6656 ,
0.53248 , 0.425984 , 0.3407872 , 0.27262976, 0.21810381])
Thanks,
Josef
>
> 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
>>
> _______________________________________________
> NumPy-Discussion mailing list
> NumPy-Discussion@scipy.org
> http://mail.scipy.org/mailman/listinfo/numpy-discussion
>
More information about the NumPy-Discussion
mailing list