1"""ARMA process and estimation with scipy.signal.lfilter 2 3Notes 4----- 5* written without textbook, works but not sure about everything 6 briefly checked and it looks to be standard least squares, see below 7 8* theoretical autocorrelation function of general ARMA 9 Done, relatively easy to guess solution, time consuming to get 10 theoretical test cases, example file contains explicit formulas for 11 acovf of MA(1), MA(2) and ARMA(1,1) 12 13Properties: 14Judge, ... (1985): The Theory and Practise of Econometrics 15 16Author: josefpktd 17License: BSD 18""" 19from statsmodels.compat.pandas import Appender 20 21import numpy as np 22from scipy import linalg, optimize, signal 23 24from statsmodels.tools.docstring import Docstring, remove_parameters 25from statsmodels.tools.validation import array_like 26 27__all__ = [ 28 "arma_acf", 29 "arma_acovf", 30 "arma_generate_sample", 31 "arma_impulse_response", 32 "arma2ar", 33 "arma2ma", 34 "deconvolve", 35 "lpol2index", 36 "index2lpol", 37] 38 39 40NONSTATIONARY_ERROR = """\ 41The model's autoregressive parameters (ar) indicate that the process 42 is non-stationary. arma_acovf can only be used with stationary processes. 43""" 44 45 46def arma_generate_sample( 47 ar, ma, nsample, scale=1, distrvs=None, axis=0, burnin=0 48): 49 """ 50 Simulate data from an ARMA. 51 52 Parameters 53 ---------- 54 ar : array_like 55 The coefficient for autoregressive lag polynomial, including zero lag. 56 ma : array_like 57 The coefficient for moving-average lag polynomial, including zero lag. 58 nsample : int or tuple of ints 59 If nsample is an integer, then this creates a 1d timeseries of 60 length size. If nsample is a tuple, creates a len(nsample) 61 dimensional time series where time is indexed along the input 62 variable ``axis``. All series are unless ``distrvs`` generates 63 dependent data. 64 scale : float 65 The standard deviation of noise. 66 distrvs : function, random number generator 67 A function that generates the random numbers, and takes ``size`` 68 as argument. The default is np.random.standard_normal. 69 axis : int 70 See nsample for details. 71 burnin : int 72 Number of observation at the beginning of the sample to drop. 73 Used to reduce dependence on initial values. 74 75 Returns 76 ------- 77 ndarray 78 Random sample(s) from an ARMA process. 79 80 Notes 81 ----- 82 As mentioned above, both the AR and MA components should include the 83 coefficient on the zero-lag. This is typically 1. Further, due to the 84 conventions used in signal processing used in signal.lfilter vs. 85 conventions in statistics for ARMA processes, the AR parameters should 86 have the opposite sign of what you might expect. See the examples below. 87 88 Examples 89 -------- 90 >>> import numpy as np 91 >>> np.random.seed(12345) 92 >>> arparams = np.array([.75, -.25]) 93 >>> maparams = np.array([.65, .35]) 94 >>> ar = np.r_[1, -arparams] # add zero-lag and negate 95 >>> ma = np.r_[1, maparams] # add zero-lag 96 >>> y = sm.tsa.arma_generate_sample(ar, ma, 250) 97 >>> model = sm.tsa.ARIMA(y, (2, 0, 2), trend='n').fit(disp=0) 98 >>> model.params 99 array([ 0.79044189, -0.23140636, 0.70072904, 0.40608028]) 100 """ 101 distrvs = np.random.standard_normal if distrvs is None else distrvs 102 if np.ndim(nsample) == 0: 103 nsample = [nsample] 104 if burnin: 105 # handle burin time for nd arrays 106 # maybe there is a better trick in scipy.fft code 107 newsize = list(nsample) 108 newsize[axis] += burnin 109 newsize = tuple(newsize) 110 fslice = [slice(None)] * len(newsize) 111 fslice[axis] = slice(burnin, None, None) 112 fslice = tuple(fslice) 113 else: 114 newsize = tuple(nsample) 115 fslice = tuple([slice(None)] * np.ndim(newsize)) 116 eta = scale * distrvs(size=newsize) 117 return signal.lfilter(ma, ar, eta, axis=axis)[fslice] 118 119 120def arma_acovf(ar, ma, nobs=10, sigma2=1, dtype=None): 121 """ 122 Theoretical autocovariances of stationary ARMA processes 123 124 Parameters 125 ---------- 126 ar : array_like, 1d 127 The coefficients for autoregressive lag polynomial, including zero lag. 128 ma : array_like, 1d 129 The coefficients for moving-average lag polynomial, including zero lag. 130 nobs : int 131 The number of terms (lags plus zero lag) to include in returned acovf. 132 sigma2 : float 133 Variance of the innovation term. 134 135 Returns 136 ------- 137 ndarray 138 The autocovariance of ARMA process given by ar, ma. 139 140 See Also 141 -------- 142 arma_acf : Autocorrelation function for ARMA processes. 143 acovf : Sample autocovariance estimation. 144 145 References 146 ---------- 147 .. [*] Brockwell, Peter J., and Richard A. Davis. 2009. Time Series: 148 Theory and Methods. 2nd ed. 1991. New York, NY: Springer. 149 """ 150 if dtype is None: 151 dtype = np.common_type(np.array(ar), np.array(ma), np.array(sigma2)) 152 153 p = len(ar) - 1 154 q = len(ma) - 1 155 m = max(p, q) + 1 156 157 if sigma2.real < 0: 158 raise ValueError("Must have positive innovation variance.") 159 160 # Short-circuit for trivial corner-case 161 if p == q == 0: 162 out = np.zeros(nobs, dtype=dtype) 163 out[0] = sigma2 164 return out 165 elif p > 0 and np.max(np.abs(np.roots(ar))) >= 1: 166 raise ValueError(NONSTATIONARY_ERROR) 167 168 # Get the moving average representation coefficients that we need 169 ma_coeffs = arma2ma(ar, ma, lags=m) 170 171 # Solve for the first m autocovariances via the linear system 172 # described by (BD, eq. 3.3.8) 173 A = np.zeros((m, m), dtype=dtype) 174 b = np.zeros((m, 1), dtype=dtype) 175 # We need a zero-right-padded version of ar params 176 tmp_ar = np.zeros(m, dtype=dtype) 177 tmp_ar[: p + 1] = ar 178 for k in range(m): 179 A[k, : (k + 1)] = tmp_ar[: (k + 1)][::-1] 180 A[k, 1 : m - k] += tmp_ar[(k + 1) : m] 181 b[k] = sigma2 * np.dot(ma[k : q + 1], ma_coeffs[: max((q + 1 - k), 0)]) 182 acovf = np.zeros(max(nobs, m), dtype=dtype) 183 try: 184 acovf[:m] = np.linalg.solve(A, b)[:, 0] 185 except np.linalg.LinAlgError: 186 raise ValueError(NONSTATIONARY_ERROR) 187 188 # Iteratively apply (BD, eq. 3.3.9) to solve for remaining autocovariances 189 if nobs > m: 190 zi = signal.lfiltic([1], ar, acovf[:m:][::-1]) 191 acovf[m:] = signal.lfilter( 192 [1], ar, np.zeros(nobs - m, dtype=dtype), zi=zi 193 )[0] 194 195 return acovf[:nobs] 196 197 198def arma_acf(ar, ma, lags=10): 199 """ 200 Theoretical autocorrelation function of an ARMA process. 201 202 Parameters 203 ---------- 204 ar : array_like 205 Coefficients for autoregressive lag polynomial, including zero lag. 206 ma : array_like 207 Coefficients for moving-average lag polynomial, including zero lag. 208 lags : int 209 The number of terms (lags plus zero lag) to include in returned acf. 210 211 Returns 212 ------- 213 ndarray 214 The autocorrelations of ARMA process given by ar and ma. 215 216 See Also 217 -------- 218 arma_acovf : Autocovariances from ARMA processes. 219 acf : Sample autocorrelation function estimation. 220 acovf : Sample autocovariance function estimation. 221 """ 222 acovf = arma_acovf(ar, ma, lags) 223 return acovf / acovf[0] 224 225 226def arma_pacf(ar, ma, lags=10): 227 """ 228 Theoretical partial autocorrelation function of an ARMA process. 229 230 Parameters 231 ---------- 232 ar : array_like, 1d 233 The coefficients for autoregressive lag polynomial, including zero lag. 234 ma : array_like, 1d 235 The coefficients for moving-average lag polynomial, including zero lag. 236 lags : int 237 The number of terms (lags plus zero lag) to include in returned pacf. 238 239 Returns 240 ------- 241 ndarrray 242 The partial autocorrelation of ARMA process given by ar and ma. 243 244 Notes 245 ----- 246 Solves yule-walker equation for each lag order up to nobs lags. 247 248 not tested/checked yet 249 """ 250 # TODO: Should use rank 1 inverse update 251 apacf = np.zeros(lags) 252 acov = arma_acf(ar, ma, lags=lags + 1) 253 254 apacf[0] = 1.0 255 for k in range(2, lags + 1): 256 r = acov[:k] 257 apacf[k - 1] = linalg.solve(linalg.toeplitz(r[:-1]), r[1:])[-1] 258 return apacf 259 260 261def arma_periodogram(ar, ma, worN=None, whole=0): 262 """ 263 Periodogram for ARMA process given by lag-polynomials ar and ma. 264 265 Parameters 266 ---------- 267 ar : array_like 268 The autoregressive lag-polynomial with leading 1 and lhs sign. 269 ma : array_like 270 The moving average lag-polynomial with leading 1. 271 worN : {None, int}, optional 272 An option for scipy.signal.freqz (read "w or N"). 273 If None, then compute at 512 frequencies around the unit circle. 274 If a single integer, the compute at that many frequencies. 275 Otherwise, compute the response at frequencies given in worN. 276 whole : {0,1}, optional 277 An options for scipy.signal.freqz/ 278 Normally, frequencies are computed from 0 to pi (upper-half of 279 unit-circle. If whole is non-zero compute frequencies from 0 to 2*pi. 280 281 Returns 282 ------- 283 w : ndarray 284 The frequencies. 285 sd : ndarray 286 The periodogram, also known as the spectral density. 287 288 Notes 289 ----- 290 Normalization ? 291 292 This uses signal.freqz, which does not use fft. There is a fft version 293 somewhere. 294 """ 295 w, h = signal.freqz(ma, ar, worN=worN, whole=whole) 296 sd = np.abs(h) ** 2 / np.sqrt(2 * np.pi) 297 if np.any(np.isnan(h)): 298 # this happens with unit root or seasonal unit root' 299 import warnings 300 301 warnings.warn( 302 "Warning: nan in frequency response h, maybe a unit " "root", 303 RuntimeWarning, 304 ) 305 return w, sd 306 307 308def arma_impulse_response(ar, ma, leads=100): 309 """ 310 Compute the impulse response function (MA representation) for ARMA process. 311 312 Parameters 313 ---------- 314 ar : array_like, 1d 315 The auto regressive lag polynomial. 316 ma : array_like, 1d 317 The moving average lag polynomial. 318 leads : int 319 The number of observations to calculate. 320 321 Returns 322 ------- 323 ndarray 324 The impulse response function with nobs elements. 325 326 Notes 327 ----- 328 This is the same as finding the MA representation of an ARMA(p,q). 329 By reversing the role of ar and ma in the function arguments, the 330 returned result is the AR representation of an ARMA(p,q), i.e 331 332 ma_representation = arma_impulse_response(ar, ma, leads=100) 333 ar_representation = arma_impulse_response(ma, ar, leads=100) 334 335 Fully tested against matlab 336 337 Examples 338 -------- 339 AR(1) 340 341 >>> arma_impulse_response([1.0, -0.8], [1.], leads=10) 342 array([ 1. , 0.8 , 0.64 , 0.512 , 0.4096 , 343 0.32768 , 0.262144 , 0.2097152 , 0.16777216, 0.13421773]) 344 345 this is the same as 346 347 >>> 0.8**np.arange(10) 348 array([ 1. , 0.8 , 0.64 , 0.512 , 0.4096 , 349 0.32768 , 0.262144 , 0.2097152 , 0.16777216, 0.13421773]) 350 351 MA(2) 352 353 >>> arma_impulse_response([1.0], [1., 0.5, 0.2], leads=10) 354 array([ 1. , 0.5, 0.2, 0. , 0. , 0. , 0. , 0. , 0. , 0. ]) 355 356 ARMA(1,2) 357 358 >>> arma_impulse_response([1.0, -0.8], [1., 0.5, 0.2], leads=10) 359 array([ 1. , 1.3 , 1.24 , 0.992 , 0.7936 , 360 0.63488 , 0.507904 , 0.4063232 , 0.32505856, 0.26004685]) 361 """ 362 impulse = np.zeros(leads) 363 impulse[0] = 1.0 364 return signal.lfilter(ma, ar, impulse) 365 366 367def arma2ma(ar, ma, lags=100): 368 """ 369 A finite-lag approximate MA representation of an ARMA process. 370 371 Parameters 372 ---------- 373 ar : ndarray 374 The auto regressive lag polynomial. 375 ma : ndarray 376 The moving average lag polynomial. 377 lags : int 378 The number of coefficients to calculate. 379 380 Returns 381 ------- 382 ndarray 383 The coefficients of AR lag polynomial with nobs elements. 384 385 Notes 386 ----- 387 Equivalent to ``arma_impulse_response(ma, ar, leads=100)`` 388 """ 389 return arma_impulse_response(ar, ma, leads=lags) 390 391 392def arma2ar(ar, ma, lags=100): 393 """ 394 A finite-lag AR approximation of an ARMA process. 395 396 Parameters 397 ---------- 398 ar : array_like 399 The auto regressive lag polynomial. 400 ma : array_like 401 The moving average lag polynomial. 402 lags : int 403 The number of coefficients to calculate. 404 405 Returns 406 ------- 407 ndarray 408 The coefficients of AR lag polynomial with nobs elements. 409 410 Notes 411 ----- 412 Equivalent to ``arma_impulse_response(ma, ar, leads=100)`` 413 """ 414 return arma_impulse_response(ma, ar, leads=lags) 415 416 417# moved from sandbox.tsa.try_fi 418def ar2arma(ar_des, p, q, n=20, mse="ar", start=None): 419 """ 420 Find arma approximation to ar process. 421 422 This finds the ARMA(p,q) coefficients that minimize the integrated 423 squared difference between the impulse_response functions (MA 424 representation) of the AR and the ARMA process. This does not check 425 whether the MA lag polynomial of the ARMA process is invertible, neither 426 does it check the roots of the AR lag polynomial. 427 428 Parameters 429 ---------- 430 ar_des : array_like 431 The coefficients of original AR lag polynomial, including lag zero. 432 p : int 433 The length of desired AR lag polynomials. 434 q : int 435 The length of desired MA lag polynomials. 436 n : int 437 The number of terms of the impulse_response function to include in the 438 objective function for the approximation. 439 mse : str, 'ar' 440 Not used. 441 start : ndarray 442 Initial values to use when finding the approximation. 443 444 Returns 445 ------- 446 ar_app : ndarray 447 The coefficients of the AR lag polynomials of the approximation. 448 ma_app : ndarray 449 The coefficients of the MA lag polynomials of the approximation. 450 res : tuple 451 The result of optimize.leastsq. 452 453 Notes 454 ----- 455 Extension is possible if we want to match autocovariance instead 456 of impulse response function. 457 """ 458 459 # TODO: convert MA lag polynomial, ma_app, to be invertible, by mirroring 460 # TODO: roots outside the unit interval to ones that are inside. How to do 461 # TODO: this? 462 463 # p,q = pq 464 def msear_err(arma, ar_des): 465 ar, ma = np.r_[1, arma[: p - 1]], np.r_[1, arma[p - 1 :]] 466 ar_approx = arma_impulse_response(ma, ar, n) 467 return ar_des - ar_approx # ((ar - ar_approx)**2).sum() 468 469 if start is None: 470 arma0 = np.r_[-0.9 * np.ones(p - 1), np.zeros(q - 1)] 471 else: 472 arma0 = start 473 res = optimize.leastsq(msear_err, arma0, ar_des, maxfev=5000) 474 arma_app = np.atleast_1d(res[0]) 475 ar_app = (np.r_[1, arma_app[: p - 1]],) 476 ma_app = np.r_[1, arma_app[p - 1 :]] 477 return ar_app, ma_app, res 478 479 480_arma_docs = {"ar": arma2ar.__doc__, "ma": arma2ma.__doc__} 481 482 483def lpol2index(ar): 484 """ 485 Remove zeros from lag polynomial 486 487 Parameters 488 ---------- 489 ar : array_like 490 coefficients of lag polynomial 491 492 Returns 493 ------- 494 coeffs : ndarray 495 non-zero coefficients of lag polynomial 496 index : ndarray 497 index (lags) of lag polynomial with non-zero elements 498 """ 499 ar = array_like(ar, "ar") 500 index = np.nonzero(ar)[0] 501 coeffs = ar[index] 502 return coeffs, index 503 504 505def index2lpol(coeffs, index): 506 """ 507 Expand coefficients to lag poly 508 509 Parameters 510 ---------- 511 coeffs : ndarray 512 non-zero coefficients of lag polynomial 513 index : ndarray 514 index (lags) of lag polynomial with non-zero elements 515 516 Returns 517 ------- 518 ar : array_like 519 coefficients of lag polynomial 520 """ 521 n = max(index) 522 ar = np.zeros(n + 1) 523 ar[index] = coeffs 524 return ar 525 526 527def lpol_fima(d, n=20): 528 """MA representation of fractional integration 529 530 .. math:: (1-L)^{-d} for |d|<0.5 or |d|<1 (?) 531 532 Parameters 533 ---------- 534 d : float 535 fractional power 536 n : int 537 number of terms to calculate, including lag zero 538 539 Returns 540 ------- 541 ma : ndarray 542 coefficients of lag polynomial 543 """ 544 # hide import inside function until we use this heavily 545 from scipy.special import gammaln 546 547 j = np.arange(n) 548 return np.exp(gammaln(d + j) - gammaln(j + 1) - gammaln(d)) 549 550 551# moved from sandbox.tsa.try_fi 552def lpol_fiar(d, n=20): 553 """AR representation of fractional integration 554 555 .. math:: (1-L)^{d} for |d|<0.5 or |d|<1 (?) 556 557 Parameters 558 ---------- 559 d : float 560 fractional power 561 n : int 562 number of terms to calculate, including lag zero 563 564 Returns 565 ------- 566 ar : ndarray 567 coefficients of lag polynomial 568 569 Notes: 570 first coefficient is 1, negative signs except for first term, 571 ar(L)*x_t 572 """ 573 # hide import inside function until we use this heavily 574 from scipy.special import gammaln 575 576 j = np.arange(n) 577 ar = -np.exp(gammaln(-d + j) - gammaln(j + 1) - gammaln(-d)) 578 ar[0] = 1 579 return ar 580 581 582# moved from sandbox.tsa.try_fi 583def lpol_sdiff(s): 584 """return coefficients for seasonal difference (1-L^s) 585 586 just a trivial convenience function 587 588 Parameters 589 ---------- 590 s : int 591 number of periods in season 592 593 Returns 594 ------- 595 sdiff : list, length s+1 596 """ 597 return [1] + [0] * (s - 1) + [-1] 598 599 600def deconvolve(num, den, n=None): 601 """Deconvolves divisor out of signal, division of polynomials for n terms 602 603 calculates den^{-1} * num 604 605 Parameters 606 ---------- 607 num : array_like 608 signal or lag polynomial 609 denom : array_like 610 coefficients of lag polynomial (linear filter) 611 n : None or int 612 number of terms of quotient 613 614 Returns 615 ------- 616 quot : ndarray 617 quotient or filtered series 618 rem : ndarray 619 remainder 620 621 Notes 622 ----- 623 If num is a time series, then this applies the linear filter den^{-1}. 624 If both num and den are both lag polynomials, then this calculates the 625 quotient polynomial for n terms and also returns the remainder. 626 627 This is copied from scipy.signal.signaltools and added n as optional 628 parameter. 629 """ 630 num = np.atleast_1d(num) 631 den = np.atleast_1d(den) 632 N = len(num) 633 D = len(den) 634 if D > N and n is None: 635 quot = [] 636 rem = num 637 else: 638 if n is None: 639 n = N - D + 1 640 input = np.zeros(n, float) 641 input[0] = 1 642 quot = signal.lfilter(num, den, input) 643 num_approx = signal.convolve(den, quot, mode="full") 644 if len(num) < len(num_approx): # 1d only ? 645 num = np.concatenate((num, np.zeros(len(num_approx) - len(num)))) 646 rem = num - num_approx 647 return quot, rem 648 649 650_generate_sample_doc = Docstring(arma_generate_sample.__doc__) 651_generate_sample_doc.remove_parameters(["ar", "ma"]) 652_generate_sample_doc.replace_block("Notes", []) 653_generate_sample_doc.replace_block("Examples", []) 654 655 656class ArmaProcess(object): 657 r""" 658 Theoretical properties of an ARMA process for specified lag-polynomials. 659 660 Parameters 661 ---------- 662 ar : array_like 663 Coefficient for autoregressive lag polynomial, including zero lag. 664 Must be entered using the signs from the lag polynomial representation. 665 See the notes for more information about the sign. 666 ma : array_like 667 Coefficient for moving-average lag polynomial, including zero lag. 668 nobs : int, optional 669 Length of simulated time series. Used, for example, if a sample is 670 generated. See example. 671 672 Notes 673 ----- 674 Both the AR and MA components must include the coefficient on the 675 zero-lag. In almost all cases these values should be 1. Further, due to 676 using the lag-polynomial representation, the AR parameters should 677 have the opposite sign of what one would write in the ARMA representation. 678 See the examples below. 679 680 The ARMA(p,q) process is described by 681 682 .. math:: 683 684 y_{t}=\phi_{1}y_{t-1}+\ldots+\phi_{p}y_{t-p}+\theta_{1}\epsilon_{t-1} 685 +\ldots+\theta_{q}\epsilon_{t-q}+\epsilon_{t} 686 687 and the parameterization used in this function uses the lag-polynomial 688 representation, 689 690 .. math:: 691 692 \left(1-\phi_{1}L-\ldots-\phi_{p}L^{p}\right)y_{t} = 693 \left(1+\theta_{1}L+\ldots+\theta_{q}L^{q}\right)\epsilon_{t} 694 695 Examples 696 -------- 697 ARMA(2,2) with AR coefficients 0.75 and -0.25, and MA coefficients 0.65 and 0.35 698 699 >>> import statsmodels.api as sm 700 >>> import numpy as np 701 >>> np.random.seed(12345) 702 >>> arparams = np.array([.75, -.25]) 703 >>> maparams = np.array([.65, .35]) 704 >>> ar = np.r_[1, -arparams] # add zero-lag and negate 705 >>> ma = np.r_[1, maparams] # add zero-lag 706 >>> arma_process = sm.tsa.ArmaProcess(ar, ma) 707 >>> arma_process.isstationary 708 True 709 >>> arma_process.isinvertible 710 True 711 >>> arma_process.arroots 712 array([1.5-1.32287566j, 1.5+1.32287566j]) 713 >>> y = arma_process.generate_sample(250) 714 >>> model = sm.tsa.ARIMA(y, (2, 0, 2), trend='n').fit(disp=0) 715 >>> model.params 716 array([ 0.79044189, -0.23140636, 0.70072904, 0.40608028]) 717 718 The same ARMA(2,2) Using the from_coeffs class method 719 720 >>> arma_process = sm.tsa.ArmaProcess.from_coeffs(arparams, maparams) 721 >>> arma_process.arroots 722 array([1.5-1.32287566j, 1.5+1.32287566j]) 723 """ 724 725 # TODO: Check unit root behavior 726 def __init__(self, ar=None, ma=None, nobs=100): 727 if ar is None: 728 ar = np.array([1.0]) 729 if ma is None: 730 ma = np.array([1.0]) 731 self.ar = array_like(ar, "ar") 732 self.ma = array_like(ma, "ma") 733 self.arcoefs = -self.ar[1:] 734 self.macoefs = self.ma[1:] 735 self.arpoly = np.polynomial.Polynomial(self.ar) 736 self.mapoly = np.polynomial.Polynomial(self.ma) 737 self.nobs = nobs 738 739 @classmethod 740 def from_roots(cls, maroots=None, arroots=None, nobs=100): 741 """ 742 Create ArmaProcess from AR and MA polynomial roots. 743 744 Parameters 745 ---------- 746 maroots : array_like 747 Roots for the MA polynomial 748 1 + theta_1*z + theta_2*z^2 + ..... + theta_n*z^n 749 arroots : array_like 750 Roots for the AR polynomial 751 1 - phi_1*z - phi_2*z^2 - ..... - phi_n*z^n 752 nobs : int, optional 753 Length of simulated time series. Used, for example, if a sample 754 is generated. 755 756 Returns 757 ------- 758 ArmaProcess 759 Class instance initialized with arcoefs and macoefs. 760 761 Examples 762 -------- 763 >>> arroots = [.75, -.25] 764 >>> maroots = [.65, .35] 765 >>> arma_process = sm.tsa.ArmaProcess.from_roots(arroots, maroots) 766 >>> arma_process.isstationary 767 True 768 >>> arma_process.isinvertible 769 True 770 """ 771 arpoly = np.polynomial.polynomial.Polynomial.fromroots(arroots) 772 mapoly = np.polynomial.polynomial.Polynomial.fromroots(maroots) 773 # As from_coeffs will create a polynomial with constant 1/-1,(MA/AR) 774 # we need to scale the polynomial coefficients accordingly 775 return cls(np.r_[1, -np.asarray(-1 * arpoly.coef[1:] / arpoly.coef[0])], 776 np.r_[1, np.asarray(mapoly.coef[1:] / mapoly.coef[0])], 777 nobs=nobs) 778 779 @classmethod 780 def from_coeffs(cls, arcoefs=None, macoefs=None, nobs=100): 781 """ 782 Create ArmaProcess from an ARMA representation. 783 784 Parameters 785 ---------- 786 arcoefs : array_like 787 Coefficient for autoregressive lag polynomial, not including zero 788 lag. The sign is inverted to conform to the usual time series 789 representation of an ARMA process in statistics. See the class 790 docstring for more information. 791 macoefs : array_like 792 Coefficient for moving-average lag polynomial, excluding zero lag. 793 nobs : int, optional 794 Length of simulated time series. Used, for example, if a sample 795 is generated. 796 797 Returns 798 ------- 799 ArmaProcess 800 Class instance initialized with arcoefs and macoefs. 801 802 Examples 803 -------- 804 >>> arparams = [.75, -.25] 805 >>> maparams = [.65, .35] 806 >>> arma_process = sm.tsa.ArmaProcess.from_coeffs(ar, ma) 807 >>> arma_process.isstationary 808 True 809 >>> arma_process.isinvertible 810 True 811 """ 812 arcoefs = [] if arcoefs is None else arcoefs 813 macoefs = [] if macoefs is None else macoefs 814 return cls( 815 np.r_[1, -np.asarray(arcoefs)], 816 np.r_[1, np.asarray(macoefs)], 817 nobs=nobs, 818 ) 819 820 @classmethod 821 def from_estimation(cls, model_results, nobs=None): 822 """ 823 Create an ArmaProcess from the results of an ARIMA estimation. 824 825 Parameters 826 ---------- 827 model_results : ARIMAResults instance 828 A fitted model. 829 nobs : int, optional 830 If None, nobs is taken from the results. 831 832 Returns 833 ------- 834 ArmaProcess 835 Class instance initialized from model_results. 836 837 See Also 838 -------- 839 statsmodels.tsa.arima.model.ARIMA 840 The models class used to create the ArmaProcess 841 """ 842 nobs = nobs or model_results.nobs 843 return cls( 844 model_results.polynomial_reduced_ar, 845 model_results.polynomial_reduced_ma, 846 nobs=nobs, 847 ) 848 849 def __mul__(self, oth): 850 if isinstance(oth, self.__class__): 851 ar = (self.arpoly * oth.arpoly).coef 852 ma = (self.mapoly * oth.mapoly).coef 853 else: 854 try: 855 aroth, maoth = oth 856 arpolyoth = np.polynomial.Polynomial(aroth) 857 mapolyoth = np.polynomial.Polynomial(maoth) 858 ar = (self.arpoly * arpolyoth).coef 859 ma = (self.mapoly * mapolyoth).coef 860 except: 861 raise TypeError("Other type is not a valid type") 862 return self.__class__(ar, ma, nobs=self.nobs) 863 864 def __repr__(self): 865 msg = "ArmaProcess({0}, {1}, nobs={2}) at {3}" 866 return msg.format( 867 self.ar.tolist(), self.ma.tolist(), self.nobs, hex(id(self)) 868 ) 869 870 def __str__(self): 871 return "ArmaProcess\nAR: {0}\nMA: {1}".format( 872 self.ar.tolist(), self.ma.tolist() 873 ) 874 875 @Appender(remove_parameters(arma_acovf.__doc__, ["ar", "ma", "sigma2"])) 876 def acovf(self, nobs=None): 877 nobs = nobs or self.nobs 878 return arma_acovf(self.ar, self.ma, nobs=nobs) 879 880 @Appender(remove_parameters(arma_acf.__doc__, ["ar", "ma"])) 881 def acf(self, lags=None): 882 lags = lags or self.nobs 883 return arma_acf(self.ar, self.ma, lags=lags) 884 885 @Appender(remove_parameters(arma_pacf.__doc__, ["ar", "ma"])) 886 def pacf(self, lags=None): 887 lags = lags or self.nobs 888 return arma_pacf(self.ar, self.ma, lags=lags) 889 890 @Appender( 891 remove_parameters( 892 arma_periodogram.__doc__, ["ar", "ma", "worN", "whole"] 893 ) 894 ) 895 def periodogram(self, nobs=None): 896 nobs = nobs or self.nobs 897 return arma_periodogram(self.ar, self.ma, worN=nobs) 898 899 @Appender(remove_parameters(arma_impulse_response.__doc__, ["ar", "ma"])) 900 def impulse_response(self, leads=None): 901 leads = leads or self.nobs 902 return arma_impulse_response(self.ar, self.ma, leads=leads) 903 904 @Appender(remove_parameters(arma2ma.__doc__, ["ar", "ma"])) 905 def arma2ma(self, lags=None): 906 lags = lags or self.lags 907 return arma2ma(self.ar, self.ma, lags=lags) 908 909 @Appender(remove_parameters(arma2ar.__doc__, ["ar", "ma"])) 910 def arma2ar(self, lags=None): 911 lags = lags or self.lags 912 return arma2ar(self.ar, self.ma, lags=lags) 913 914 @property 915 def arroots(self): 916 """Roots of autoregressive lag-polynomial""" 917 return self.arpoly.roots() 918 919 @property 920 def maroots(self): 921 """Roots of moving average lag-polynomial""" 922 return self.mapoly.roots() 923 924 @property 925 def isstationary(self): 926 """ 927 Arma process is stationary if AR roots are outside unit circle. 928 929 Returns 930 ------- 931 bool 932 True if autoregressive roots are outside unit circle. 933 """ 934 if np.all(np.abs(self.arroots) > 1.0): 935 return True 936 else: 937 return False 938 939 @property 940 def isinvertible(self): 941 """ 942 Arma process is invertible if MA roots are outside unit circle. 943 944 Returns 945 ------- 946 bool 947 True if moving average roots are outside unit circle. 948 """ 949 if np.all(np.abs(self.maroots) > 1): 950 return True 951 else: 952 return False 953 954 def invertroots(self, retnew=False): 955 """ 956 Make MA polynomial invertible by inverting roots inside unit circle. 957 958 Parameters 959 ---------- 960 retnew : bool 961 If False (default), then return the lag-polynomial as array. 962 If True, then return a new instance with invertible MA-polynomial. 963 964 Returns 965 ------- 966 manew : ndarray 967 A new invertible MA lag-polynomial, returned if retnew is false. 968 wasinvertible : bool 969 True if the MA lag-polynomial was already invertible, returned if 970 retnew is false. 971 armaprocess : new instance of class 972 If retnew is true, then return a new instance with invertible 973 MA-polynomial. 974 """ 975 # TODO: variable returns like this? 976 pr = self.maroots 977 mainv = self.ma 978 invertible = self.isinvertible 979 if not invertible: 980 pr[np.abs(pr) < 1] = 1.0 / pr[np.abs(pr) < 1] 981 pnew = np.polynomial.Polynomial.fromroots(pr) 982 mainv = pnew.coef / pnew.coef[0] 983 984 if retnew: 985 return self.__class__(self.ar, mainv, nobs=self.nobs) 986 else: 987 return mainv, invertible 988 989 @Appender(str(_generate_sample_doc)) 990 def generate_sample( 991 self, nsample=100, scale=1.0, distrvs=None, axis=0, burnin=0 992 ): 993 return arma_generate_sample( 994 self.ar, self.ma, nsample, scale, distrvs, axis=axis, burnin=burnin 995 ) 996