1#cython: boundscheck=False 2#cython: wraparound=False 3#cython: cdivision=False 4""" 5Innovations algorithm 6 7Author: Chad Fulton 8License: Simplified-BSD 9""" 10 11# Typical imports 12import numpy as np 13import warnings 14from statsmodels.tsa import arima_process 15from statsmodels.src.math cimport NPY_PI, dlog, zlog 16cimport numpy as cnp 17cimport cython 18 19cnp.import_array() 20 21cdef int C = 0 22 23 24 25cdef stoeplitz(int n, int offset0, int offset1, 26 cnp.float32_t [:] in_column, 27 cnp.float32_t [:, :] out_matrix): 28 """ 29 toeplitz(int n, int o0, int o1, cnp.float32_t [:] in_column, cnp.float32_t [:, :] out_matrix 30 31 Construct a Toeplitz block in a matrix in place 32 33 Parameters 34 ---------- 35 n : int 36 The number of entries of `in_column` to use. 37 offset0 : int 38 The row offset for `out_matrix` at which to begin writing the Toeplitz 39 block. 40 offset1 : int 41 The column offset for `out_matrix` at which to begin writing the 42 Toeplitz block. 43 in_column : ndarray 44 The column used to construct the Toeplitz block. 45 out_matrix : ndarray 46 The matrix in which to write a Toeplitz block 47 48 Notes 49 ----- 50 This function constructs the Toeplitz block in-place and does not return 51 any output. 52 """ 53 cdef Py_ssize_t i, j 54 55 for i in range(n): 56 for j in range(i + 1): 57 out_matrix[offset0 + i, offset1 + j] = in_column[i - j] 58 if i != j: 59 # Note: scipy by default does complex conjugate, but not 60 # necessary here since we're only dealing with covariances, 61 # which will be real (except in the case of complex-step 62 # differentiation, but in that case we do not want to apply 63 # the conjugation anyway) 64 out_matrix[offset0 + j, offset1 + i] = in_column[i - j] 65 66 67cpdef sarma_transformed_acovf_fast(cnp.float32_t [:] ar, 68 cnp.float32_t [:] ma, 69 cnp.float32_t [:] arma_acovf): 70 """ 71 arma_transformed_acovf_fast(cnp.float32_t [:] ar, cnp.float32_t [:] ma, cnp.float32_t [:] arma_acovf) 72 73 Quickly construct the autocovariance matrix for a transformed process. 74 75 Using the autocovariance function for an ARMA process, constructs the 76 autocovariances associated with the transformed process described 77 in equation (3.3.1) of _[1] in a memory efficient, and so fast, way. 78 79 Parameters 80 ---------- 81 ar : ndarray 82 Autoregressive coefficients, including the zero lag, where the sign 83 convention assumes the coefficients are part of the lag polynomial on 84 the left-hand-side of the ARMA definition (i.e. they have the opposite 85 sign from the usual econometrics convention in which the coefficients 86 are on the right-hand-side of the ARMA definition). 87 ma : ndarray 88 Moving average coefficients, including the zero lag, where the sign 89 convention assumes the coefficients are part of the lag polynomial on 90 the right-hand-side of the ARMA definition (i.e. they have the same 91 sign from the usual econometrics convention in which the coefficients 92 are on the right-hand-side of the ARMA definition). 93 arma_acovf : ndarray 94 The vector of autocovariances of the ARMA process. 95 96 Returns 97 ------- 98 acovf : ndarray 99 A matrix containing the autocovariances of the portion of the 100 transformed process with time-varying autocovariances. Its dimension 101 is equal to `min(m * 2, n)` where `m = max(len(ar) - 1, len(ma) - 1)` 102 and `n` is the length of the input array `arma_acovf`. It is important 103 to note that only the values in the first `m` columns or `m` rows are 104 valid. In particular, the entries in the block `acovf[m:, m:]` should 105 not be used in any case (and in fact will always be equal to zeros). 106 acovf2 : ndarray 107 An array containing the autocovariance function of the portion of the 108 transformed process with time-invariant autocovariances. Its dimension 109 is equal to `max(n - m, 0)` where `n` is the length of the input 110 array `arma_acovf`. 111 112 Notes 113 ----- 114 The definition of this autocovariance matrix is from _[1] equation 3.3.3. 115 116 This function assumes that the variance of the ARMA innovation term is 117 equal to one. If this is not true, then the calling program should replace 118 `arma_acovf` with `arma_acovf / sigma2`, where sigma2 is that variance. 119 120 This function is relatively fast even when `arma_acovf` is large, since 121 it only constructs the full autocovariance matrix for a generally small 122 subset of observations. The trade-off is that the output of this function 123 is somewhat more difficult to use. 124 125 References 126 ---------- 127 .. [1] Brockwell, Peter J., and Richard A. Davis. 2009. 128 Time Series: Theory and Methods. 2nd ed. 1991. 129 New York, NY: Springer. 130 """ 131 cdef Py_ssize_t nobs, p, q, m, m2, n, i, j, r, tmp_ix 132 cdef cnp.npy_intp dim1[1] 133 cdef cnp.npy_intp dim2[2] 134 cdef cnp.float32_t [:, :] acovf 135 cdef cnp.float32_t [:] acovf2 136 137 nobs = arma_acovf.shape[0] 138 p = len(ar) - 1 139 q = len(ma) - 1 140 m = max(p, q) 141 m2 = 2 * m 142 n = min(m2, nobs) 143 144 dim2[0] = m2; 145 dim2[1] = m2; 146 acovf = cnp.PyArray_ZEROS(2, dim2, cnp.NPY_FLOAT32, C) 147 dim1[0] = max(nobs - m, 0); 148 acovf2 = cnp.PyArray_ZEROS(1, dim1, cnp.NPY_FLOAT32, C) 149 150 # Set i, j = 1, ..., m (which is then done) 151 stoeplitz(m, 0, 0, arma_acovf, acovf) 152 153 # Set i = 1, ..., m; j = m + 1, ..., 2 m 154 # and i = m + 1, ..., 2 m; j = 1, ..., m 155 if nobs > m: 156 for j in range(m): 157 for i in range(m, m2): 158 acovf[i, j] = arma_acovf[i - j] 159 for r in range(1, p + 1): 160 tmp_ix = abs(r - (i - j)) 161 acovf[i, j] = acovf[i, j] - (-ar[r] * arma_acovf[tmp_ix]) 162 acovf[:m, m:m2] = acovf[m:m2, :m].T 163 164 # Set values for |i - j| <= q, min(i, j) = m + 1, and max(i, j) <= nobs 165 if nobs > m: 166 for i in range(nobs - m): 167 for r in range(q + 1 - i): 168 acovf2[i] = acovf2[i] + ma[r] * ma[r + i] 169 170 return acovf[:n, :n], acovf2 171 172 173cpdef sarma_innovations_algo_fast(int nobs, 174 cnp.float32_t [:] ar_params, 175 cnp.float32_t [:] ma_params, 176 cnp.float32_t [:, :] acovf, 177 cnp.float32_t [:] acovf2): 178 """ 179 arma_innovations_algo_fast(int nobs, cnp.float32_t [:] ar_params, cnp.float32_t [:] ma_params, cnp.float32_t [:, :] acovf, cnp.float32_t [:] acovf2) 180 181 Quickly apply innovations algorithm for an ARMA process. 182 183 Parameters 184 ---------- 185 ar_params : ndarray 186 Autoregressive parameters. 187 ma_params : ndarray 188 Moving average parameters. 189 acovf : ndarray 190 An `m * 2` x `m * 2` autocovariance matrix at least the first `m` 191 columns filled in, where `m = max(len(ar_params), ma_params)` 192 (see `arma_transformed_acovf_fast`). 193 acovf2 : ndarray 194 A `max(0, nobs - m)` length vector containing the autocovariance 195 function associated with the final `nobs - m` observations. 196 197 Returns 198 ------- 199 theta : ndarray 200 The matrix of moving average coefficients from the innovations 201 algorithm. 202 v : ndarray 203 The vector of mean squared errors. 204 205 Notes 206 ----- 207 The innovations algorithm is presented in _[1], section 2.5.2. 208 209 Details of the innovations algorithm applied to ARMA processes is given 210 in _[1] section 3.3 and in section 5.2.7. 211 212 This function is relatively fast even with a large number of observations 213 since we can exploit a number of known zeros in the theta array. 214 215 References 216 ---------- 217 .. [1] Brockwell, Peter J., and Richard A. Davis. 2009. 218 Time Series: Theory and Methods. 2nd ed. 1991. 219 New York, NY: Springer. 220 """ 221 cdef Py_ssize_t i, j, k, n, _n, m, m2, p, q, start, start2 222 cdef cnp.npy_intp dim1[1] 223 cdef cnp.npy_intp dim2[2] 224 cdef cnp.float32_t [:, :] theta 225 cdef cnp.float32_t [:] v 226 227 p = len(ar_params) 228 q = len(ma_params) 229 m = max(p, q) 230 m2 = 2 * m 231 232 dim1[0] = nobs; 233 v = cnp.PyArray_ZEROS(1, dim1, cnp.NPY_FLOAT32, C) 234 dim2[0] = nobs; 235 dim2[1] = m + 1; 236 theta = cnp.PyArray_ZEROS(2, dim2, cnp.NPY_FLOAT32, C) 237 238 if m > 0: 239 v[0] = acovf[0, 0] 240 else: 241 # (handle the corner case where p = q = 0) 242 v[0] = acovf2[0] 243 244 for n in range(nobs - 1): 245 _n = n + 1 246 247 start = 0 if n < m else n + 1 - q 248 for k in range(start, n + 1): 249 # See Brockwell and Davis, p. 100-101 250 # (here we have weak rather than strict inequality due to Python's 251 # zero indexing) 252 if n >= m and n - k >= q: 253 continue 254 255 if n + 1 < m2 and k < m: 256 theta[_n, n - k] = acovf[n + 1, k] 257 else: 258 theta[_n, n - k] = acovf2[n + 1 - k] 259 260 start2 = 0 if n < m else n - m 261 for j in range(start2, k): 262 if n - j < m + 1: 263 theta[_n, n - k] = theta[_n, n - k] - theta[k - 1 + 1, k - j - 1] * theta[_n, n - j] * v[j] 264 theta[_n, n - k] = theta[_n, n - k] / v[k] 265 266 if n + 1 < m: 267 v[n + 1] = acovf[n + 1, n + 1] 268 else: 269 v[n + 1] = acovf2[0] 270 start = max(0, n - (m + 1) + 2) 271 for i in range(start, n + 1): 272 v[n + 1] = v[n + 1] - theta[_n, n - i]**2 * v[i] 273 274 return theta, v 275 276 277cpdef sarma_innovations_filter(cnp.float32_t [:] endog, 278 cnp.float32_t [:] ar_params, 279 cnp.float32_t [:] ma_params, 280 cnp.float32_t [:, :] theta): 281 """ 282 arma_innovations_filter(cnp.float32_t [:] endog, cnp.float32_t [:] ar_params, cnp.float32_t [:] ma_params, cnp.float32_t [:, :] theta): 283 284 Innovations filter for an ARMA process. 285 286 Parameters 287 ---------- 288 endog : ndarray 289 The observed time-series process. 290 ar_params : ndarray 291 Autoregressive parameters. 292 ma_params : ndarray 293 Moving average parameters. 294 theta : ndarray 295 The matrix of moving average coefficients from the innovations 296 algorithm (see `arma_innovations_algo` or `arma_innovations_algo_fast`) 297 298 Returns 299 ------- 300 u : ndarray 301 The vector of innovations: the one-step-ahead prediction errors from 302 applying the innovations algorithm. 303 304 Notes 305 ----- 306 The innovations algorithm is presented in _[1], section 2.5.2. 307 308 Details of the innovations algorithm applied to ARMA processes is given 309 in _[1] section 3.3 and in section 5.2.7. 310 311 References 312 ---------- 313 .. [1] Brockwell, Peter J., and Richard A. Davis. 2009. 314 Time Series: Theory and Methods. 2nd ed. 1991. 315 New York, NY: Springer. 316 """ 317 cdef Py_ssize_t nobs, i, k, j, m, p, q 318 cdef cnp.npy_intp dim1[1] 319 cdef cnp.float32_t [:] u 320 cdef cnp.float32_t hat 321 322 p = len(ar_params) 323 q = len(ma_params) 324 m = max(p, q) 325 326 nobs = theta.shape[0] 327 k = theta.shape[1] 328 329 dim1[0] = nobs; 330 u = cnp.PyArray_ZEROS(1, dim1, cnp.NPY_FLOAT32, C) 331 u[0] = endog[0] 332 333 for i in range(1, nobs): 334 hat = 0 335 if i < m: 336 for j in range(i): 337 hat = hat + theta[i, j] * u[i - j - 1] 338 else: 339 for j in range(p): 340 hat = hat + ar_params[j] * endog[i - j - 1] 341 for j in range(q): 342 hat = hat + theta[i, j] * u[i - j - 1] 343 u[i] = endog[i] - hat 344 345 return u 346 347 348cpdef sarma_innovations(cnp.float32_t [:] endog, 349 cnp.float32_t [:] ar_params, 350 cnp.float32_t [:] ma_params, 351 cnp.float32_t sigma2): 352 """ 353 arma_innovations(cnp.float32_t [:] endog, cnp.float32_t [:] ar_params, cnp.float32_t [:] ma_params): 354 355 Compute innovations and variances based on an ARMA process. 356 357 Parameters 358 ---------- 359 endog : ndarray 360 The observed time-series process. 361 ar_params : ndarray 362 Autoregressive parameters. 363 ma_params : ndarray 364 Moving average parameters. 365 sigma2 : ndarray 366 The ARMA innovation variance. 367 368 Returns 369 ------- 370 u : ndarray 371 The vector of innovations: the one-step-ahead prediction errors from 372 applying the innovations algorithm. 373 v : ndarray 374 The vector of innovation variances. 375 376 Notes 377 ----- 378 The innovations algorithm is presented in _[1], section 2.5.2. 379 380 References 381 ---------- 382 .. [1] Brockwell, Peter J., and Richard A. Davis. 2009. 383 Time Series: Theory and Methods. 2nd ed. 1991. 384 New York, NY: Springer. 385 """ 386 cdef Py_ssize_t nobs = len(endog), i 387 cdef cnp.float32_t const 388 cdef cnp.float32_t [:] ar, ma, arma_acovf, llf_obs, acovf2, u, v 389 cdef cnp.float32_t [:, :] acovf 390 cdef cnp.npy_intp dim1[1] 391 392 dim1[0] = len(ar_params) + 1; 393 ar = cnp.PyArray_ZEROS(1, dim1, cnp.NPY_FLOAT32, C) 394 dim1[0] = len(ma_params) + 1; 395 ma = cnp.PyArray_ZEROS(1, dim1, cnp.NPY_FLOAT32, C) 396 397 ar[0] = 1 398 for i in range(1, len(ar_params) + 1): 399 ar[i] = -1 * ar_params[i - 1] 400 ma[0] = 1 401 ma[1:] = ma_params 402 403 arma_acovf = arima_process.arma_acovf(ar, ma, nobs, sigma2, dtype=np.float32) / sigma2 404 acovf, acovf2 = sarma_transformed_acovf_fast(ar, ma, arma_acovf) 405 theta, v = sarma_innovations_algo_fast(nobs, ar_params, ma_params, acovf, acovf2) 406 u = sarma_innovations_filter(endog, ar_params, ma_params, theta) 407 408 return u, v 409 410 411cpdef sarma_loglikeobs_fast(cnp.float32_t [:] endog, 412 cnp.float32_t [:] ar_params, 413 cnp.float32_t [:] ma_params, 414 cnp.float32_t sigma2): 415 """ 416 sarma_loglikeobs_fast(cnp.float32_t [:] endog, cnp.float32_t [:] ar_params, cnp.float32_t [:] ma_params, cnp.float32_t sigma2) 417 418 Quickly calculate the loglikelihood of each observation for an ARMA process 419 420 Parameters 421 ---------- 422 endog : ndarray 423 The observed time-series process. 424 ar_params : ndarray 425 Autoregressive parameters. 426 ma_params : ndarray 427 Moving average parameters. 428 sigma2 : ndarray 429 The ARMA innovation variance. 430 431 Returns 432 ------- 433 loglike : ndarray of float 434 Array of loglikelihood values for each observation. 435 436 Notes 437 ----- 438 Details related to computing the loglikelihood associated with an ARMA 439 process using the innovations algorithm are given in _[1] section 5.2. 440 441 References 442 ---------- 443 .. [1] Brockwell, Peter J., and Richard A. Davis. 2009. 444 Time Series: Theory and Methods. 2nd ed. 1991. 445 New York, NY: Springer. 446 """ 447 448 cdef Py_ssize_t nobs = len(endog), i 449 cdef cnp.float32_t const 450 cdef cnp.float32_t [:] llf_obs, u, v 451 cdef cnp.npy_intp dim1[1] 452 453 u, v = sarma_innovations(endog, ar_params, ma_params, sigma2) 454 455 dim1[0] = nobs; 456 llf_obs = cnp.PyArray_ZEROS(1, dim1, cnp.NPY_FLOAT32, C) 457 458 const = dlog(2*NPY_PI) 459 for i in range(nobs): 460 llf_obs[i] = -0.5 * u[i]**2 / (sigma2 * v[i]) - 0.5 * (const + dlog(sigma2 * v[i])) 461 462 return np.array(llf_obs, dtype=np.float32) 463 464 465cdef dtoeplitz(int n, int offset0, int offset1, 466 cnp.float64_t [:] in_column, 467 cnp.float64_t [:, :] out_matrix): 468 """ 469 toeplitz(int n, int o0, int o1, cnp.float64_t [:] in_column, cnp.float64_t [:, :] out_matrix 470 471 Construct a Toeplitz block in a matrix in place 472 473 Parameters 474 ---------- 475 n : int 476 The number of entries of `in_column` to use. 477 offset0 : int 478 The row offset for `out_matrix` at which to begin writing the Toeplitz 479 block. 480 offset1 : int 481 The column offset for `out_matrix` at which to begin writing the 482 Toeplitz block. 483 in_column : ndarray 484 The column used to construct the Toeplitz block. 485 out_matrix : ndarray 486 The matrix in which to write a Toeplitz block 487 488 Notes 489 ----- 490 This function constructs the Toeplitz block in-place and does not return 491 any output. 492 """ 493 cdef Py_ssize_t i, j 494 495 for i in range(n): 496 for j in range(i + 1): 497 out_matrix[offset0 + i, offset1 + j] = in_column[i - j] 498 if i != j: 499 # Note: scipy by default does complex conjugate, but not 500 # necessary here since we're only dealing with covariances, 501 # which will be real (except in the case of complex-step 502 # differentiation, but in that case we do not want to apply 503 # the conjugation anyway) 504 out_matrix[offset0 + j, offset1 + i] = in_column[i - j] 505 506 507cpdef darma_transformed_acovf_fast(cnp.float64_t [:] ar, 508 cnp.float64_t [:] ma, 509 cnp.float64_t [:] arma_acovf): 510 """ 511 arma_transformed_acovf_fast(cnp.float64_t [:] ar, cnp.float64_t [:] ma, cnp.float64_t [:] arma_acovf) 512 513 Quickly construct the autocovariance matrix for a transformed process. 514 515 Using the autocovariance function for an ARMA process, constructs the 516 autocovariances associated with the transformed process described 517 in equation (3.3.1) of _[1] in a memory efficient, and so fast, way. 518 519 Parameters 520 ---------- 521 ar : ndarray 522 Autoregressive coefficients, including the zero lag, where the sign 523 convention assumes the coefficients are part of the lag polynomial on 524 the left-hand-side of the ARMA definition (i.e. they have the opposite 525 sign from the usual econometrics convention in which the coefficients 526 are on the right-hand-side of the ARMA definition). 527 ma : ndarray 528 Moving average coefficients, including the zero lag, where the sign 529 convention assumes the coefficients are part of the lag polynomial on 530 the right-hand-side of the ARMA definition (i.e. they have the same 531 sign from the usual econometrics convention in which the coefficients 532 are on the right-hand-side of the ARMA definition). 533 arma_acovf : ndarray 534 The vector of autocovariances of the ARMA process. 535 536 Returns 537 ------- 538 acovf : ndarray 539 A matrix containing the autocovariances of the portion of the 540 transformed process with time-varying autocovariances. Its dimension 541 is equal to `min(m * 2, n)` where `m = max(len(ar) - 1, len(ma) - 1)` 542 and `n` is the length of the input array `arma_acovf`. It is important 543 to note that only the values in the first `m` columns or `m` rows are 544 valid. In particular, the entries in the block `acovf[m:, m:]` should 545 not be used in any case (and in fact will always be equal to zeros). 546 acovf2 : ndarray 547 An array containing the autocovariance function of the portion of the 548 transformed process with time-invariant autocovariances. Its dimension 549 is equal to `max(n - m, 0)` where `n` is the length of the input 550 array `arma_acovf`. 551 552 Notes 553 ----- 554 The definition of this autocovariance matrix is from _[1] equation 3.3.3. 555 556 This function assumes that the variance of the ARMA innovation term is 557 equal to one. If this is not true, then the calling program should replace 558 `arma_acovf` with `arma_acovf / sigma2`, where sigma2 is that variance. 559 560 This function is relatively fast even when `arma_acovf` is large, since 561 it only constructs the full autocovariance matrix for a generally small 562 subset of observations. The trade-off is that the output of this function 563 is somewhat more difficult to use. 564 565 References 566 ---------- 567 .. [1] Brockwell, Peter J., and Richard A. Davis. 2009. 568 Time Series: Theory and Methods. 2nd ed. 1991. 569 New York, NY: Springer. 570 """ 571 cdef Py_ssize_t nobs, p, q, m, m2, n, i, j, r, tmp_ix 572 cdef cnp.npy_intp dim1[1] 573 cdef cnp.npy_intp dim2[2] 574 cdef cnp.float64_t [:, :] acovf 575 cdef cnp.float64_t [:] acovf2 576 577 nobs = arma_acovf.shape[0] 578 p = len(ar) - 1 579 q = len(ma) - 1 580 m = max(p, q) 581 m2 = 2 * m 582 n = min(m2, nobs) 583 584 dim2[0] = m2; 585 dim2[1] = m2; 586 acovf = cnp.PyArray_ZEROS(2, dim2, cnp.NPY_FLOAT64, C) 587 dim1[0] = max(nobs - m, 0); 588 acovf2 = cnp.PyArray_ZEROS(1, dim1, cnp.NPY_FLOAT64, C) 589 590 # Set i, j = 1, ..., m (which is then done) 591 dtoeplitz(m, 0, 0, arma_acovf, acovf) 592 593 # Set i = 1, ..., m; j = m + 1, ..., 2 m 594 # and i = m + 1, ..., 2 m; j = 1, ..., m 595 if nobs > m: 596 for j in range(m): 597 for i in range(m, m2): 598 acovf[i, j] = arma_acovf[i - j] 599 for r in range(1, p + 1): 600 tmp_ix = abs(r - (i - j)) 601 acovf[i, j] = acovf[i, j] - (-ar[r] * arma_acovf[tmp_ix]) 602 acovf[:m, m:m2] = acovf[m:m2, :m].T 603 604 # Set values for |i - j| <= q, min(i, j) = m + 1, and max(i, j) <= nobs 605 if nobs > m: 606 for i in range(nobs - m): 607 for r in range(q + 1 - i): 608 acovf2[i] = acovf2[i] + ma[r] * ma[r + i] 609 610 return acovf[:n, :n], acovf2 611 612 613cpdef darma_innovations_algo_fast(int nobs, 614 cnp.float64_t [:] ar_params, 615 cnp.float64_t [:] ma_params, 616 cnp.float64_t [:, :] acovf, 617 cnp.float64_t [:] acovf2): 618 """ 619 arma_innovations_algo_fast(int nobs, cnp.float64_t [:] ar_params, cnp.float64_t [:] ma_params, cnp.float64_t [:, :] acovf, cnp.float64_t [:] acovf2) 620 621 Quickly apply innovations algorithm for an ARMA process. 622 623 Parameters 624 ---------- 625 ar_params : ndarray 626 Autoregressive parameters. 627 ma_params : ndarray 628 Moving average parameters. 629 acovf : ndarray 630 An `m * 2` x `m * 2` autocovariance matrix at least the first `m` 631 columns filled in, where `m = max(len(ar_params), ma_params)` 632 (see `arma_transformed_acovf_fast`). 633 acovf2 : ndarray 634 A `max(0, nobs - m)` length vector containing the autocovariance 635 function associated with the final `nobs - m` observations. 636 637 Returns 638 ------- 639 theta : ndarray 640 The matrix of moving average coefficients from the innovations 641 algorithm. 642 v : ndarray 643 The vector of mean squared errors. 644 645 Notes 646 ----- 647 The innovations algorithm is presented in _[1], section 2.5.2. 648 649 Details of the innovations algorithm applied to ARMA processes is given 650 in _[1] section 3.3 and in section 5.2.7. 651 652 This function is relatively fast even with a large number of observations 653 since we can exploit a number of known zeros in the theta array. 654 655 References 656 ---------- 657 .. [1] Brockwell, Peter J., and Richard A. Davis. 2009. 658 Time Series: Theory and Methods. 2nd ed. 1991. 659 New York, NY: Springer. 660 """ 661 cdef Py_ssize_t i, j, k, n, _n, m, m2, p, q, start, start2 662 cdef cnp.npy_intp dim1[1] 663 cdef cnp.npy_intp dim2[2] 664 cdef cnp.float64_t [:, :] theta 665 cdef cnp.float64_t [:] v 666 667 p = len(ar_params) 668 q = len(ma_params) 669 m = max(p, q) 670 m2 = 2 * m 671 672 dim1[0] = nobs; 673 v = cnp.PyArray_ZEROS(1, dim1, cnp.NPY_FLOAT64, C) 674 dim2[0] = nobs; 675 dim2[1] = m + 1; 676 theta = cnp.PyArray_ZEROS(2, dim2, cnp.NPY_FLOAT64, C) 677 678 if m > 0: 679 v[0] = acovf[0, 0] 680 else: 681 # (handle the corner case where p = q = 0) 682 v[0] = acovf2[0] 683 684 for n in range(nobs - 1): 685 _n = n + 1 686 687 start = 0 if n < m else n + 1 - q 688 for k in range(start, n + 1): 689 # See Brockwell and Davis, p. 100-101 690 # (here we have weak rather than strict inequality due to Python's 691 # zero indexing) 692 if n >= m and n - k >= q: 693 continue 694 695 if n + 1 < m2 and k < m: 696 theta[_n, n - k] = acovf[n + 1, k] 697 else: 698 theta[_n, n - k] = acovf2[n + 1 - k] 699 700 start2 = 0 if n < m else n - m 701 for j in range(start2, k): 702 if n - j < m + 1: 703 theta[_n, n - k] = theta[_n, n - k] - theta[k - 1 + 1, k - j - 1] * theta[_n, n - j] * v[j] 704 theta[_n, n - k] = theta[_n, n - k] / v[k] 705 706 if n + 1 < m: 707 v[n + 1] = acovf[n + 1, n + 1] 708 else: 709 v[n + 1] = acovf2[0] 710 start = max(0, n - (m + 1) + 2) 711 for i in range(start, n + 1): 712 v[n + 1] = v[n + 1] - theta[_n, n - i]**2 * v[i] 713 714 return theta, v 715 716 717cpdef darma_innovations_filter(cnp.float64_t [:] endog, 718 cnp.float64_t [:] ar_params, 719 cnp.float64_t [:] ma_params, 720 cnp.float64_t [:, :] theta): 721 """ 722 arma_innovations_filter(cnp.float64_t [:] endog, cnp.float64_t [:] ar_params, cnp.float64_t [:] ma_params, cnp.float64_t [:, :] theta): 723 724 Innovations filter for an ARMA process. 725 726 Parameters 727 ---------- 728 endog : ndarray 729 The observed time-series process. 730 ar_params : ndarray 731 Autoregressive parameters. 732 ma_params : ndarray 733 Moving average parameters. 734 theta : ndarray 735 The matrix of moving average coefficients from the innovations 736 algorithm (see `arma_innovations_algo` or `arma_innovations_algo_fast`) 737 738 Returns 739 ------- 740 u : ndarray 741 The vector of innovations: the one-step-ahead prediction errors from 742 applying the innovations algorithm. 743 744 Notes 745 ----- 746 The innovations algorithm is presented in _[1], section 2.5.2. 747 748 Details of the innovations algorithm applied to ARMA processes is given 749 in _[1] section 3.3 and in section 5.2.7. 750 751 References 752 ---------- 753 .. [1] Brockwell, Peter J., and Richard A. Davis. 2009. 754 Time Series: Theory and Methods. 2nd ed. 1991. 755 New York, NY: Springer. 756 """ 757 cdef Py_ssize_t nobs, i, k, j, m, p, q 758 cdef cnp.npy_intp dim1[1] 759 cdef cnp.float64_t [:] u 760 cdef cnp.float64_t hat 761 762 p = len(ar_params) 763 q = len(ma_params) 764 m = max(p, q) 765 766 nobs = theta.shape[0] 767 k = theta.shape[1] 768 769 dim1[0] = nobs; 770 u = cnp.PyArray_ZEROS(1, dim1, cnp.NPY_FLOAT64, C) 771 u[0] = endog[0] 772 773 for i in range(1, nobs): 774 hat = 0 775 if i < m: 776 for j in range(i): 777 hat = hat + theta[i, j] * u[i - j - 1] 778 else: 779 for j in range(p): 780 hat = hat + ar_params[j] * endog[i - j - 1] 781 for j in range(q): 782 hat = hat + theta[i, j] * u[i - j - 1] 783 u[i] = endog[i] - hat 784 785 return u 786 787 788cpdef darma_innovations(cnp.float64_t [:] endog, 789 cnp.float64_t [:] ar_params, 790 cnp.float64_t [:] ma_params, 791 cnp.float64_t sigma2): 792 """ 793 arma_innovations(cnp.float64_t [:] endog, cnp.float64_t [:] ar_params, cnp.float64_t [:] ma_params): 794 795 Compute innovations and variances based on an ARMA process. 796 797 Parameters 798 ---------- 799 endog : ndarray 800 The observed time-series process. 801 ar_params : ndarray 802 Autoregressive parameters. 803 ma_params : ndarray 804 Moving average parameters. 805 sigma2 : ndarray 806 The ARMA innovation variance. 807 808 Returns 809 ------- 810 u : ndarray 811 The vector of innovations: the one-step-ahead prediction errors from 812 applying the innovations algorithm. 813 v : ndarray 814 The vector of innovation variances. 815 816 Notes 817 ----- 818 The innovations algorithm is presented in _[1], section 2.5.2. 819 820 References 821 ---------- 822 .. [1] Brockwell, Peter J., and Richard A. Davis. 2009. 823 Time Series: Theory and Methods. 2nd ed. 1991. 824 New York, NY: Springer. 825 """ 826 cdef Py_ssize_t nobs = len(endog), i 827 cdef cnp.float64_t const 828 cdef cnp.float64_t [:] ar, ma, arma_acovf, llf_obs, acovf2, u, v 829 cdef cnp.float64_t [:, :] acovf 830 cdef cnp.npy_intp dim1[1] 831 832 dim1[0] = len(ar_params) + 1; 833 ar = cnp.PyArray_ZEROS(1, dim1, cnp.NPY_FLOAT64, C) 834 dim1[0] = len(ma_params) + 1; 835 ma = cnp.PyArray_ZEROS(1, dim1, cnp.NPY_FLOAT64, C) 836 837 ar[0] = 1 838 for i in range(1, len(ar_params) + 1): 839 ar[i] = -1 * ar_params[i - 1] 840 ma[0] = 1 841 ma[1:] = ma_params 842 843 arma_acovf = arima_process.arma_acovf(ar, ma, nobs, sigma2, dtype=float) / sigma2 844 acovf, acovf2 = darma_transformed_acovf_fast(ar, ma, arma_acovf) 845 theta, v = darma_innovations_algo_fast(nobs, ar_params, ma_params, acovf, acovf2) 846 u = darma_innovations_filter(endog, ar_params, ma_params, theta) 847 848 return u, v 849 850 851cpdef darma_loglikeobs_fast(cnp.float64_t [:] endog, 852 cnp.float64_t [:] ar_params, 853 cnp.float64_t [:] ma_params, 854 cnp.float64_t sigma2): 855 """ 856 darma_loglikeobs_fast(cnp.float64_t [:] endog, cnp.float64_t [:] ar_params, cnp.float64_t [:] ma_params, cnp.float64_t sigma2) 857 858 Quickly calculate the loglikelihood of each observation for an ARMA process 859 860 Parameters 861 ---------- 862 endog : ndarray 863 The observed time-series process. 864 ar_params : ndarray 865 Autoregressive parameters. 866 ma_params : ndarray 867 Moving average parameters. 868 sigma2 : ndarray 869 The ARMA innovation variance. 870 871 Returns 872 ------- 873 loglike : ndarray of float 874 Array of loglikelihood values for each observation. 875 876 Notes 877 ----- 878 Details related to computing the loglikelihood associated with an ARMA 879 process using the innovations algorithm are given in _[1] section 5.2. 880 881 References 882 ---------- 883 .. [1] Brockwell, Peter J., and Richard A. Davis. 2009. 884 Time Series: Theory and Methods. 2nd ed. 1991. 885 New York, NY: Springer. 886 """ 887 888 cdef Py_ssize_t nobs = len(endog), i 889 cdef cnp.float64_t const 890 cdef cnp.float64_t [:] llf_obs, u, v 891 cdef cnp.npy_intp dim1[1] 892 893 u, v = darma_innovations(endog, ar_params, ma_params, sigma2) 894 895 dim1[0] = nobs; 896 llf_obs = cnp.PyArray_ZEROS(1, dim1, cnp.NPY_FLOAT64, C) 897 898 const = dlog(2*NPY_PI) 899 for i in range(nobs): 900 llf_obs[i] = -0.5 * u[i]**2 / (sigma2 * v[i]) - 0.5 * (const + dlog(sigma2 * v[i])) 901 902 return np.array(llf_obs, dtype=float) 903 904 905cdef ctoeplitz(int n, int offset0, int offset1, 906 cnp.complex64_t [:] in_column, 907 cnp.complex64_t [:, :] out_matrix): 908 """ 909 toeplitz(int n, int o0, int o1, cnp.complex64_t [:] in_column, cnp.complex64_t [:, :] out_matrix 910 911 Construct a Toeplitz block in a matrix in place 912 913 Parameters 914 ---------- 915 n : int 916 The number of entries of `in_column` to use. 917 offset0 : int 918 The row offset for `out_matrix` at which to begin writing the Toeplitz 919 block. 920 offset1 : int 921 The column offset for `out_matrix` at which to begin writing the 922 Toeplitz block. 923 in_column : ndarray 924 The column used to construct the Toeplitz block. 925 out_matrix : ndarray 926 The matrix in which to write a Toeplitz block 927 928 Notes 929 ----- 930 This function constructs the Toeplitz block in-place and does not return 931 any output. 932 """ 933 cdef Py_ssize_t i, j 934 935 for i in range(n): 936 for j in range(i + 1): 937 out_matrix[offset0 + i, offset1 + j] = in_column[i - j] 938 if i != j: 939 # Note: scipy by default does complex conjugate, but not 940 # necessary here since we're only dealing with covariances, 941 # which will be real (except in the case of complex-step 942 # differentiation, but in that case we do not want to apply 943 # the conjugation anyway) 944 out_matrix[offset0 + j, offset1 + i] = in_column[i - j] 945 946 947cpdef carma_transformed_acovf_fast(cnp.complex64_t [:] ar, 948 cnp.complex64_t [:] ma, 949 cnp.complex64_t [:] arma_acovf): 950 """ 951 arma_transformed_acovf_fast(cnp.complex64_t [:] ar, cnp.complex64_t [:] ma, cnp.complex64_t [:] arma_acovf) 952 953 Quickly construct the autocovariance matrix for a transformed process. 954 955 Using the autocovariance function for an ARMA process, constructs the 956 autocovariances associated with the transformed process described 957 in equation (3.3.1) of _[1] in a memory efficient, and so fast, way. 958 959 Parameters 960 ---------- 961 ar : ndarray 962 Autoregressive coefficients, including the zero lag, where the sign 963 convention assumes the coefficients are part of the lag polynomial on 964 the left-hand-side of the ARMA definition (i.e. they have the opposite 965 sign from the usual econometrics convention in which the coefficients 966 are on the right-hand-side of the ARMA definition). 967 ma : ndarray 968 Moving average coefficients, including the zero lag, where the sign 969 convention assumes the coefficients are part of the lag polynomial on 970 the right-hand-side of the ARMA definition (i.e. they have the same 971 sign from the usual econometrics convention in which the coefficients 972 are on the right-hand-side of the ARMA definition). 973 arma_acovf : ndarray 974 The vector of autocovariances of the ARMA process. 975 976 Returns 977 ------- 978 acovf : ndarray 979 A matrix containing the autocovariances of the portion of the 980 transformed process with time-varying autocovariances. Its dimension 981 is equal to `min(m * 2, n)` where `m = max(len(ar) - 1, len(ma) - 1)` 982 and `n` is the length of the input array `arma_acovf`. It is important 983 to note that only the values in the first `m` columns or `m` rows are 984 valid. In particular, the entries in the block `acovf[m:, m:]` should 985 not be used in any case (and in fact will always be equal to zeros). 986 acovf2 : ndarray 987 An array containing the autocovariance function of the portion of the 988 transformed process with time-invariant autocovariances. Its dimension 989 is equal to `max(n - m, 0)` where `n` is the length of the input 990 array `arma_acovf`. 991 992 Notes 993 ----- 994 The definition of this autocovariance matrix is from _[1] equation 3.3.3. 995 996 This function assumes that the variance of the ARMA innovation term is 997 equal to one. If this is not true, then the calling program should replace 998 `arma_acovf` with `arma_acovf / sigma2`, where sigma2 is that variance. 999 1000 This function is relatively fast even when `arma_acovf` is large, since 1001 it only constructs the full autocovariance matrix for a generally small 1002 subset of observations. The trade-off is that the output of this function 1003 is somewhat more difficult to use. 1004 1005 References 1006 ---------- 1007 .. [1] Brockwell, Peter J., and Richard A. Davis. 2009. 1008 Time Series: Theory and Methods. 2nd ed. 1991. 1009 New York, NY: Springer. 1010 """ 1011 cdef Py_ssize_t nobs, p, q, m, m2, n, i, j, r, tmp_ix 1012 cdef cnp.npy_intp dim1[1] 1013 cdef cnp.npy_intp dim2[2] 1014 cdef cnp.complex64_t [:, :] acovf 1015 cdef cnp.complex64_t [:] acovf2 1016 1017 nobs = arma_acovf.shape[0] 1018 p = len(ar) - 1 1019 q = len(ma) - 1 1020 m = max(p, q) 1021 m2 = 2 * m 1022 n = min(m2, nobs) 1023 1024 dim2[0] = m2; 1025 dim2[1] = m2; 1026 acovf = cnp.PyArray_ZEROS(2, dim2, cnp.NPY_COMPLEX64, C) 1027 dim1[0] = max(nobs - m, 0); 1028 acovf2 = cnp.PyArray_ZEROS(1, dim1, cnp.NPY_COMPLEX64, C) 1029 1030 # Set i, j = 1, ..., m (which is then done) 1031 ctoeplitz(m, 0, 0, arma_acovf, acovf) 1032 1033 # Set i = 1, ..., m; j = m + 1, ..., 2 m 1034 # and i = m + 1, ..., 2 m; j = 1, ..., m 1035 if nobs > m: 1036 for j in range(m): 1037 for i in range(m, m2): 1038 acovf[i, j] = arma_acovf[i - j] 1039 for r in range(1, p + 1): 1040 tmp_ix = abs(r - (i - j)) 1041 acovf[i, j] = acovf[i, j] - (-ar[r] * arma_acovf[tmp_ix]) 1042 acovf[:m, m:m2] = acovf[m:m2, :m].T 1043 1044 # Set values for |i - j| <= q, min(i, j) = m + 1, and max(i, j) <= nobs 1045 if nobs > m: 1046 for i in range(nobs - m): 1047 for r in range(q + 1 - i): 1048 acovf2[i] = acovf2[i] + ma[r] * ma[r + i] 1049 1050 return acovf[:n, :n], acovf2 1051 1052 1053cpdef carma_innovations_algo_fast(int nobs, 1054 cnp.complex64_t [:] ar_params, 1055 cnp.complex64_t [:] ma_params, 1056 cnp.complex64_t [:, :] acovf, 1057 cnp.complex64_t [:] acovf2): 1058 """ 1059 arma_innovations_algo_fast(int nobs, cnp.complex64_t [:] ar_params, cnp.complex64_t [:] ma_params, cnp.complex64_t [:, :] acovf, cnp.complex64_t [:] acovf2) 1060 1061 Quickly apply innovations algorithm for an ARMA process. 1062 1063 Parameters 1064 ---------- 1065 ar_params : ndarray 1066 Autoregressive parameters. 1067 ma_params : ndarray 1068 Moving average parameters. 1069 acovf : ndarray 1070 An `m * 2` x `m * 2` autocovariance matrix at least the first `m` 1071 columns filled in, where `m = max(len(ar_params), ma_params)` 1072 (see `arma_transformed_acovf_fast`). 1073 acovf2 : ndarray 1074 A `max(0, nobs - m)` length vector containing the autocovariance 1075 function associated with the final `nobs - m` observations. 1076 1077 Returns 1078 ------- 1079 theta : ndarray 1080 The matrix of moving average coefficients from the innovations 1081 algorithm. 1082 v : ndarray 1083 The vector of mean squared errors. 1084 1085 Notes 1086 ----- 1087 The innovations algorithm is presented in _[1], section 2.5.2. 1088 1089 Details of the innovations algorithm applied to ARMA processes is given 1090 in _[1] section 3.3 and in section 5.2.7. 1091 1092 This function is relatively fast even with a large number of observations 1093 since we can exploit a number of known zeros in the theta array. 1094 1095 References 1096 ---------- 1097 .. [1] Brockwell, Peter J., and Richard A. Davis. 2009. 1098 Time Series: Theory and Methods. 2nd ed. 1991. 1099 New York, NY: Springer. 1100 """ 1101 cdef Py_ssize_t i, j, k, n, _n, m, m2, p, q, start, start2 1102 cdef cnp.npy_intp dim1[1] 1103 cdef cnp.npy_intp dim2[2] 1104 cdef cnp.complex64_t [:, :] theta 1105 cdef cnp.complex64_t [:] v 1106 1107 p = len(ar_params) 1108 q = len(ma_params) 1109 m = max(p, q) 1110 m2 = 2 * m 1111 1112 dim1[0] = nobs; 1113 v = cnp.PyArray_ZEROS(1, dim1, cnp.NPY_COMPLEX64, C) 1114 dim2[0] = nobs; 1115 dim2[1] = m + 1; 1116 theta = cnp.PyArray_ZEROS(2, dim2, cnp.NPY_COMPLEX64, C) 1117 1118 if m > 0: 1119 v[0] = acovf[0, 0] 1120 else: 1121 # (handle the corner case where p = q = 0) 1122 v[0] = acovf2[0] 1123 1124 for n in range(nobs - 1): 1125 _n = n + 1 1126 1127 start = 0 if n < m else n + 1 - q 1128 for k in range(start, n + 1): 1129 # See Brockwell and Davis, p. 100-101 1130 # (here we have weak rather than strict inequality due to Python's 1131 # zero indexing) 1132 if n >= m and n - k >= q: 1133 continue 1134 1135 if n + 1 < m2 and k < m: 1136 theta[_n, n - k] = acovf[n + 1, k] 1137 else: 1138 theta[_n, n - k] = acovf2[n + 1 - k] 1139 1140 start2 = 0 if n < m else n - m 1141 for j in range(start2, k): 1142 if n - j < m + 1: 1143 theta[_n, n - k] = theta[_n, n - k] - theta[k - 1 + 1, k - j - 1] * theta[_n, n - j] * v[j] 1144 theta[_n, n - k] = theta[_n, n - k] / v[k] 1145 1146 if n + 1 < m: 1147 v[n + 1] = acovf[n + 1, n + 1] 1148 else: 1149 v[n + 1] = acovf2[0] 1150 start = max(0, n - (m + 1) + 2) 1151 for i in range(start, n + 1): 1152 v[n + 1] = v[n + 1] - theta[_n, n - i]**2 * v[i] 1153 1154 return theta, v 1155 1156 1157cpdef carma_innovations_filter(cnp.complex64_t [:] endog, 1158 cnp.complex64_t [:] ar_params, 1159 cnp.complex64_t [:] ma_params, 1160 cnp.complex64_t [:, :] theta): 1161 """ 1162 arma_innovations_filter(cnp.complex64_t [:] endog, cnp.complex64_t [:] ar_params, cnp.complex64_t [:] ma_params, cnp.complex64_t [:, :] theta): 1163 1164 Innovations filter for an ARMA process. 1165 1166 Parameters 1167 ---------- 1168 endog : ndarray 1169 The observed time-series process. 1170 ar_params : ndarray 1171 Autoregressive parameters. 1172 ma_params : ndarray 1173 Moving average parameters. 1174 theta : ndarray 1175 The matrix of moving average coefficients from the innovations 1176 algorithm (see `arma_innovations_algo` or `arma_innovations_algo_fast`) 1177 1178 Returns 1179 ------- 1180 u : ndarray 1181 The vector of innovations: the one-step-ahead prediction errors from 1182 applying the innovations algorithm. 1183 1184 Notes 1185 ----- 1186 The innovations algorithm is presented in _[1], section 2.5.2. 1187 1188 Details of the innovations algorithm applied to ARMA processes is given 1189 in _[1] section 3.3 and in section 5.2.7. 1190 1191 References 1192 ---------- 1193 .. [1] Brockwell, Peter J., and Richard A. Davis. 2009. 1194 Time Series: Theory and Methods. 2nd ed. 1991. 1195 New York, NY: Springer. 1196 """ 1197 cdef Py_ssize_t nobs, i, k, j, m, p, q 1198 cdef cnp.npy_intp dim1[1] 1199 cdef cnp.complex64_t [:] u 1200 cdef cnp.complex64_t hat 1201 1202 p = len(ar_params) 1203 q = len(ma_params) 1204 m = max(p, q) 1205 1206 nobs = theta.shape[0] 1207 k = theta.shape[1] 1208 1209 dim1[0] = nobs; 1210 u = cnp.PyArray_ZEROS(1, dim1, cnp.NPY_COMPLEX64, C) 1211 u[0] = endog[0] 1212 1213 for i in range(1, nobs): 1214 hat = 0 1215 if i < m: 1216 for j in range(i): 1217 hat = hat + theta[i, j] * u[i - j - 1] 1218 else: 1219 for j in range(p): 1220 hat = hat + ar_params[j] * endog[i - j - 1] 1221 for j in range(q): 1222 hat = hat + theta[i, j] * u[i - j - 1] 1223 u[i] = endog[i] - hat 1224 1225 return u 1226 1227 1228cpdef carma_innovations(cnp.complex64_t [:] endog, 1229 cnp.complex64_t [:] ar_params, 1230 cnp.complex64_t [:] ma_params, 1231 cnp.complex64_t sigma2): 1232 """ 1233 arma_innovations(cnp.complex64_t [:] endog, cnp.complex64_t [:] ar_params, cnp.complex64_t [:] ma_params): 1234 1235 Compute innovations and variances based on an ARMA process. 1236 1237 Parameters 1238 ---------- 1239 endog : ndarray 1240 The observed time-series process. 1241 ar_params : ndarray 1242 Autoregressive parameters. 1243 ma_params : ndarray 1244 Moving average parameters. 1245 sigma2 : ndarray 1246 The ARMA innovation variance. 1247 1248 Returns 1249 ------- 1250 u : ndarray 1251 The vector of innovations: the one-step-ahead prediction errors from 1252 applying the innovations algorithm. 1253 v : ndarray 1254 The vector of innovation variances. 1255 1256 Notes 1257 ----- 1258 The innovations algorithm is presented in _[1], section 2.5.2. 1259 1260 References 1261 ---------- 1262 .. [1] Brockwell, Peter J., and Richard A. Davis. 2009. 1263 Time Series: Theory and Methods. 2nd ed. 1991. 1264 New York, NY: Springer. 1265 """ 1266 cdef Py_ssize_t nobs = len(endog), i 1267 cdef cnp.complex64_t const 1268 cdef cnp.complex64_t [:] ar, ma, arma_acovf, llf_obs, acovf2, u, v 1269 cdef cnp.complex64_t [:, :] acovf 1270 cdef cnp.npy_intp dim1[1] 1271 1272 dim1[0] = len(ar_params) + 1; 1273 ar = cnp.PyArray_ZEROS(1, dim1, cnp.NPY_COMPLEX64, C) 1274 dim1[0] = len(ma_params) + 1; 1275 ma = cnp.PyArray_ZEROS(1, dim1, cnp.NPY_COMPLEX64, C) 1276 1277 ar[0] = 1 1278 for i in range(1, len(ar_params) + 1): 1279 ar[i] = -1 * ar_params[i - 1] 1280 ma[0] = 1 1281 ma[1:] = ma_params 1282 1283 arma_acovf = arima_process.arma_acovf(ar, ma, nobs, sigma2, dtype=np.complex64) / sigma2 1284 acovf, acovf2 = carma_transformed_acovf_fast(ar, ma, arma_acovf) 1285 theta, v = carma_innovations_algo_fast(nobs, ar_params, ma_params, acovf, acovf2) 1286 u = carma_innovations_filter(endog, ar_params, ma_params, theta) 1287 1288 return u, v 1289 1290 1291cpdef carma_loglikeobs_fast(cnp.complex64_t [:] endog, 1292 cnp.complex64_t [:] ar_params, 1293 cnp.complex64_t [:] ma_params, 1294 cnp.complex64_t sigma2): 1295 """ 1296 carma_loglikeobs_fast(cnp.complex64_t [:] endog, cnp.complex64_t [:] ar_params, cnp.complex64_t [:] ma_params, cnp.complex64_t sigma2) 1297 1298 Quickly calculate the loglikelihood of each observation for an ARMA process 1299 1300 Parameters 1301 ---------- 1302 endog : ndarray 1303 The observed time-series process. 1304 ar_params : ndarray 1305 Autoregressive parameters. 1306 ma_params : ndarray 1307 Moving average parameters. 1308 sigma2 : ndarray 1309 The ARMA innovation variance. 1310 1311 Returns 1312 ------- 1313 loglike : ndarray of float 1314 Array of loglikelihood values for each observation. 1315 1316 Notes 1317 ----- 1318 Details related to computing the loglikelihood associated with an ARMA 1319 process using the innovations algorithm are given in _[1] section 5.2. 1320 1321 References 1322 ---------- 1323 .. [1] Brockwell, Peter J., and Richard A. Davis. 2009. 1324 Time Series: Theory and Methods. 2nd ed. 1991. 1325 New York, NY: Springer. 1326 """ 1327 1328 cdef Py_ssize_t nobs = len(endog), i 1329 cdef cnp.complex64_t const 1330 cdef cnp.complex64_t [:] llf_obs, u, v 1331 cdef cnp.npy_intp dim1[1] 1332 1333 u, v = carma_innovations(endog, ar_params, ma_params, sigma2) 1334 1335 dim1[0] = nobs; 1336 llf_obs = cnp.PyArray_ZEROS(1, dim1, cnp.NPY_COMPLEX64, C) 1337 1338 const = zlog(2*NPY_PI) 1339 for i in range(nobs): 1340 llf_obs[i] = -0.5 * u[i]**2 / (sigma2 * v[i]) - 0.5 * (const + zlog(sigma2 * v[i])) 1341 1342 return np.array(llf_obs, dtype=np.complex64) 1343 1344 1345cdef ztoeplitz(int n, int offset0, int offset1, 1346 cnp.complex128_t [:] in_column, 1347 cnp.complex128_t [:, :] out_matrix): 1348 """ 1349 toeplitz(int n, int o0, int o1, cnp.complex128_t [:] in_column, cnp.complex128_t [:, :] out_matrix 1350 1351 Construct a Toeplitz block in a matrix in place 1352 1353 Parameters 1354 ---------- 1355 n : int 1356 The number of entries of `in_column` to use. 1357 offset0 : int 1358 The row offset for `out_matrix` at which to begin writing the Toeplitz 1359 block. 1360 offset1 : int 1361 The column offset for `out_matrix` at which to begin writing the 1362 Toeplitz block. 1363 in_column : ndarray 1364 The column used to construct the Toeplitz block. 1365 out_matrix : ndarray 1366 The matrix in which to write a Toeplitz block 1367 1368 Notes 1369 ----- 1370 This function constructs the Toeplitz block in-place and does not return 1371 any output. 1372 """ 1373 cdef Py_ssize_t i, j 1374 1375 for i in range(n): 1376 for j in range(i + 1): 1377 out_matrix[offset0 + i, offset1 + j] = in_column[i - j] 1378 if i != j: 1379 # Note: scipy by default does complex conjugate, but not 1380 # necessary here since we're only dealing with covariances, 1381 # which will be real (except in the case of complex-step 1382 # differentiation, but in that case we do not want to apply 1383 # the conjugation anyway) 1384 out_matrix[offset0 + j, offset1 + i] = in_column[i - j] 1385 1386 1387cpdef zarma_transformed_acovf_fast(cnp.complex128_t [:] ar, 1388 cnp.complex128_t [:] ma, 1389 cnp.complex128_t [:] arma_acovf): 1390 """ 1391 arma_transformed_acovf_fast(cnp.complex128_t [:] ar, cnp.complex128_t [:] ma, cnp.complex128_t [:] arma_acovf) 1392 1393 Quickly construct the autocovariance matrix for a transformed process. 1394 1395 Using the autocovariance function for an ARMA process, constructs the 1396 autocovariances associated with the transformed process described 1397 in equation (3.3.1) of _[1] in a memory efficient, and so fast, way. 1398 1399 Parameters 1400 ---------- 1401 ar : ndarray 1402 Autoregressive coefficients, including the zero lag, where the sign 1403 convention assumes the coefficients are part of the lag polynomial on 1404 the left-hand-side of the ARMA definition (i.e. they have the opposite 1405 sign from the usual econometrics convention in which the coefficients 1406 are on the right-hand-side of the ARMA definition). 1407 ma : ndarray 1408 Moving average coefficients, including the zero lag, where the sign 1409 convention assumes the coefficients are part of the lag polynomial on 1410 the right-hand-side of the ARMA definition (i.e. they have the same 1411 sign from the usual econometrics convention in which the coefficients 1412 are on the right-hand-side of the ARMA definition). 1413 arma_acovf : ndarray 1414 The vector of autocovariances of the ARMA process. 1415 1416 Returns 1417 ------- 1418 acovf : ndarray 1419 A matrix containing the autocovariances of the portion of the 1420 transformed process with time-varying autocovariances. Its dimension 1421 is equal to `min(m * 2, n)` where `m = max(len(ar) - 1, len(ma) - 1)` 1422 and `n` is the length of the input array `arma_acovf`. It is important 1423 to note that only the values in the first `m` columns or `m` rows are 1424 valid. In particular, the entries in the block `acovf[m:, m:]` should 1425 not be used in any case (and in fact will always be equal to zeros). 1426 acovf2 : ndarray 1427 An array containing the autocovariance function of the portion of the 1428 transformed process with time-invariant autocovariances. Its dimension 1429 is equal to `max(n - m, 0)` where `n` is the length of the input 1430 array `arma_acovf`. 1431 1432 Notes 1433 ----- 1434 The definition of this autocovariance matrix is from _[1] equation 3.3.3. 1435 1436 This function assumes that the variance of the ARMA innovation term is 1437 equal to one. If this is not true, then the calling program should replace 1438 `arma_acovf` with `arma_acovf / sigma2`, where sigma2 is that variance. 1439 1440 This function is relatively fast even when `arma_acovf` is large, since 1441 it only constructs the full autocovariance matrix for a generally small 1442 subset of observations. The trade-off is that the output of this function 1443 is somewhat more difficult to use. 1444 1445 References 1446 ---------- 1447 .. [1] Brockwell, Peter J., and Richard A. Davis. 2009. 1448 Time Series: Theory and Methods. 2nd ed. 1991. 1449 New York, NY: Springer. 1450 """ 1451 cdef Py_ssize_t nobs, p, q, m, m2, n, i, j, r, tmp_ix 1452 cdef cnp.npy_intp dim1[1] 1453 cdef cnp.npy_intp dim2[2] 1454 cdef cnp.complex128_t [:, :] acovf 1455 cdef cnp.complex128_t [:] acovf2 1456 1457 nobs = arma_acovf.shape[0] 1458 p = len(ar) - 1 1459 q = len(ma) - 1 1460 m = max(p, q) 1461 m2 = 2 * m 1462 n = min(m2, nobs) 1463 1464 dim2[0] = m2; 1465 dim2[1] = m2; 1466 acovf = cnp.PyArray_ZEROS(2, dim2, cnp.NPY_COMPLEX128, C) 1467 dim1[0] = max(nobs - m, 0); 1468 acovf2 = cnp.PyArray_ZEROS(1, dim1, cnp.NPY_COMPLEX128, C) 1469 1470 # Set i, j = 1, ..., m (which is then done) 1471 ztoeplitz(m, 0, 0, arma_acovf, acovf) 1472 1473 # Set i = 1, ..., m; j = m + 1, ..., 2 m 1474 # and i = m + 1, ..., 2 m; j = 1, ..., m 1475 if nobs > m: 1476 for j in range(m): 1477 for i in range(m, m2): 1478 acovf[i, j] = arma_acovf[i - j] 1479 for r in range(1, p + 1): 1480 tmp_ix = abs(r - (i - j)) 1481 acovf[i, j] = acovf[i, j] - (-ar[r] * arma_acovf[tmp_ix]) 1482 acovf[:m, m:m2] = acovf[m:m2, :m].T 1483 1484 # Set values for |i - j| <= q, min(i, j) = m + 1, and max(i, j) <= nobs 1485 if nobs > m: 1486 for i in range(nobs - m): 1487 for r in range(q + 1 - i): 1488 acovf2[i] = acovf2[i] + ma[r] * ma[r + i] 1489 1490 return acovf[:n, :n], acovf2 1491 1492 1493cpdef zarma_innovations_algo_fast(int nobs, 1494 cnp.complex128_t [:] ar_params, 1495 cnp.complex128_t [:] ma_params, 1496 cnp.complex128_t [:, :] acovf, 1497 cnp.complex128_t [:] acovf2): 1498 """ 1499 arma_innovations_algo_fast(int nobs, cnp.complex128_t [:] ar_params, cnp.complex128_t [:] ma_params, cnp.complex128_t [:, :] acovf, cnp.complex128_t [:] acovf2) 1500 1501 Quickly apply innovations algorithm for an ARMA process. 1502 1503 Parameters 1504 ---------- 1505 ar_params : ndarray 1506 Autoregressive parameters. 1507 ma_params : ndarray 1508 Moving average parameters. 1509 acovf : ndarray 1510 An `m * 2` x `m * 2` autocovariance matrix at least the first `m` 1511 columns filled in, where `m = max(len(ar_params), ma_params)` 1512 (see `arma_transformed_acovf_fast`). 1513 acovf2 : ndarray 1514 A `max(0, nobs - m)` length vector containing the autocovariance 1515 function associated with the final `nobs - m` observations. 1516 1517 Returns 1518 ------- 1519 theta : ndarray 1520 The matrix of moving average coefficients from the innovations 1521 algorithm. 1522 v : ndarray 1523 The vector of mean squared errors. 1524 1525 Notes 1526 ----- 1527 The innovations algorithm is presented in _[1], section 2.5.2. 1528 1529 Details of the innovations algorithm applied to ARMA processes is given 1530 in _[1] section 3.3 and in section 5.2.7. 1531 1532 This function is relatively fast even with a large number of observations 1533 since we can exploit a number of known zeros in the theta array. 1534 1535 References 1536 ---------- 1537 .. [1] Brockwell, Peter J., and Richard A. Davis. 2009. 1538 Time Series: Theory and Methods. 2nd ed. 1991. 1539 New York, NY: Springer. 1540 """ 1541 cdef Py_ssize_t i, j, k, n, _n, m, m2, p, q, start, start2 1542 cdef cnp.npy_intp dim1[1] 1543 cdef cnp.npy_intp dim2[2] 1544 cdef cnp.complex128_t [:, :] theta 1545 cdef cnp.complex128_t [:] v 1546 1547 p = len(ar_params) 1548 q = len(ma_params) 1549 m = max(p, q) 1550 m2 = 2 * m 1551 1552 dim1[0] = nobs; 1553 v = cnp.PyArray_ZEROS(1, dim1, cnp.NPY_COMPLEX128, C) 1554 dim2[0] = nobs; 1555 dim2[1] = m + 1; 1556 theta = cnp.PyArray_ZEROS(2, dim2, cnp.NPY_COMPLEX128, C) 1557 1558 if m > 0: 1559 v[0] = acovf[0, 0] 1560 else: 1561 # (handle the corner case where p = q = 0) 1562 v[0] = acovf2[0] 1563 1564 for n in range(nobs - 1): 1565 _n = n + 1 1566 1567 start = 0 if n < m else n + 1 - q 1568 for k in range(start, n + 1): 1569 # See Brockwell and Davis, p. 100-101 1570 # (here we have weak rather than strict inequality due to Python's 1571 # zero indexing) 1572 if n >= m and n - k >= q: 1573 continue 1574 1575 if n + 1 < m2 and k < m: 1576 theta[_n, n - k] = acovf[n + 1, k] 1577 else: 1578 theta[_n, n - k] = acovf2[n + 1 - k] 1579 1580 start2 = 0 if n < m else n - m 1581 for j in range(start2, k): 1582 if n - j < m + 1: 1583 theta[_n, n - k] = theta[_n, n - k] - theta[k - 1 + 1, k - j - 1] * theta[_n, n - j] * v[j] 1584 theta[_n, n - k] = theta[_n, n - k] / v[k] 1585 1586 if n + 1 < m: 1587 v[n + 1] = acovf[n + 1, n + 1] 1588 else: 1589 v[n + 1] = acovf2[0] 1590 start = max(0, n - (m + 1) + 2) 1591 for i in range(start, n + 1): 1592 v[n + 1] = v[n + 1] - theta[_n, n - i]**2 * v[i] 1593 1594 return theta, v 1595 1596 1597cpdef zarma_innovations_filter(cnp.complex128_t [:] endog, 1598 cnp.complex128_t [:] ar_params, 1599 cnp.complex128_t [:] ma_params, 1600 cnp.complex128_t [:, :] theta): 1601 """ 1602 arma_innovations_filter(cnp.complex128_t [:] endog, cnp.complex128_t [:] ar_params, cnp.complex128_t [:] ma_params, cnp.complex128_t [:, :] theta): 1603 1604 Innovations filter for an ARMA process. 1605 1606 Parameters 1607 ---------- 1608 endog : ndarray 1609 The observed time-series process. 1610 ar_params : ndarray 1611 Autoregressive parameters. 1612 ma_params : ndarray 1613 Moving average parameters. 1614 theta : ndarray 1615 The matrix of moving average coefficients from the innovations 1616 algorithm (see `arma_innovations_algo` or `arma_innovations_algo_fast`) 1617 1618 Returns 1619 ------- 1620 u : ndarray 1621 The vector of innovations: the one-step-ahead prediction errors from 1622 applying the innovations algorithm. 1623 1624 Notes 1625 ----- 1626 The innovations algorithm is presented in _[1], section 2.5.2. 1627 1628 Details of the innovations algorithm applied to ARMA processes is given 1629 in _[1] section 3.3 and in section 5.2.7. 1630 1631 References 1632 ---------- 1633 .. [1] Brockwell, Peter J., and Richard A. Davis. 2009. 1634 Time Series: Theory and Methods. 2nd ed. 1991. 1635 New York, NY: Springer. 1636 """ 1637 cdef Py_ssize_t nobs, i, k, j, m, p, q 1638 cdef cnp.npy_intp dim1[1] 1639 cdef cnp.complex128_t [:] u 1640 cdef cnp.complex128_t hat 1641 1642 p = len(ar_params) 1643 q = len(ma_params) 1644 m = max(p, q) 1645 1646 nobs = theta.shape[0] 1647 k = theta.shape[1] 1648 1649 dim1[0] = nobs; 1650 u = cnp.PyArray_ZEROS(1, dim1, cnp.NPY_COMPLEX128, C) 1651 u[0] = endog[0] 1652 1653 for i in range(1, nobs): 1654 hat = 0 1655 if i < m: 1656 for j in range(i): 1657 hat = hat + theta[i, j] * u[i - j - 1] 1658 else: 1659 for j in range(p): 1660 hat = hat + ar_params[j] * endog[i - j - 1] 1661 for j in range(q): 1662 hat = hat + theta[i, j] * u[i - j - 1] 1663 u[i] = endog[i] - hat 1664 1665 return u 1666 1667 1668cpdef zarma_innovations(cnp.complex128_t [:] endog, 1669 cnp.complex128_t [:] ar_params, 1670 cnp.complex128_t [:] ma_params, 1671 cnp.complex128_t sigma2): 1672 """ 1673 arma_innovations(cnp.complex128_t [:] endog, cnp.complex128_t [:] ar_params, cnp.complex128_t [:] ma_params): 1674 1675 Compute innovations and variances based on an ARMA process. 1676 1677 Parameters 1678 ---------- 1679 endog : ndarray 1680 The observed time-series process. 1681 ar_params : ndarray 1682 Autoregressive parameters. 1683 ma_params : ndarray 1684 Moving average parameters. 1685 sigma2 : ndarray 1686 The ARMA innovation variance. 1687 1688 Returns 1689 ------- 1690 u : ndarray 1691 The vector of innovations: the one-step-ahead prediction errors from 1692 applying the innovations algorithm. 1693 v : ndarray 1694 The vector of innovation variances. 1695 1696 Notes 1697 ----- 1698 The innovations algorithm is presented in _[1], section 2.5.2. 1699 1700 References 1701 ---------- 1702 .. [1] Brockwell, Peter J., and Richard A. Davis. 2009. 1703 Time Series: Theory and Methods. 2nd ed. 1991. 1704 New York, NY: Springer. 1705 """ 1706 cdef Py_ssize_t nobs = len(endog), i 1707 cdef cnp.complex128_t const 1708 cdef cnp.complex128_t [:] ar, ma, arma_acovf, llf_obs, acovf2, u, v 1709 cdef cnp.complex128_t [:, :] acovf 1710 cdef cnp.npy_intp dim1[1] 1711 1712 dim1[0] = len(ar_params) + 1; 1713 ar = cnp.PyArray_ZEROS(1, dim1, cnp.NPY_COMPLEX128, C) 1714 dim1[0] = len(ma_params) + 1; 1715 ma = cnp.PyArray_ZEROS(1, dim1, cnp.NPY_COMPLEX128, C) 1716 1717 ar[0] = 1 1718 for i in range(1, len(ar_params) + 1): 1719 ar[i] = -1 * ar_params[i - 1] 1720 ma[0] = 1 1721 ma[1:] = ma_params 1722 1723 arma_acovf = arima_process.arma_acovf(ar, ma, nobs, sigma2, dtype=complex) / sigma2 1724 acovf, acovf2 = zarma_transformed_acovf_fast(ar, ma, arma_acovf) 1725 theta, v = zarma_innovations_algo_fast(nobs, ar_params, ma_params, acovf, acovf2) 1726 u = zarma_innovations_filter(endog, ar_params, ma_params, theta) 1727 1728 return u, v 1729 1730 1731cpdef zarma_loglikeobs_fast(cnp.complex128_t [:] endog, 1732 cnp.complex128_t [:] ar_params, 1733 cnp.complex128_t [:] ma_params, 1734 cnp.complex128_t sigma2): 1735 """ 1736 zarma_loglikeobs_fast(cnp.complex128_t [:] endog, cnp.complex128_t [:] ar_params, cnp.complex128_t [:] ma_params, cnp.complex128_t sigma2) 1737 1738 Quickly calculate the loglikelihood of each observation for an ARMA process 1739 1740 Parameters 1741 ---------- 1742 endog : ndarray 1743 The observed time-series process. 1744 ar_params : ndarray 1745 Autoregressive parameters. 1746 ma_params : ndarray 1747 Moving average parameters. 1748 sigma2 : ndarray 1749 The ARMA innovation variance. 1750 1751 Returns 1752 ------- 1753 loglike : ndarray of float 1754 Array of loglikelihood values for each observation. 1755 1756 Notes 1757 ----- 1758 Details related to computing the loglikelihood associated with an ARMA 1759 process using the innovations algorithm are given in _[1] section 5.2. 1760 1761 References 1762 ---------- 1763 .. [1] Brockwell, Peter J., and Richard A. Davis. 2009. 1764 Time Series: Theory and Methods. 2nd ed. 1991. 1765 New York, NY: Springer. 1766 """ 1767 1768 cdef Py_ssize_t nobs = len(endog), i 1769 cdef cnp.complex128_t const 1770 cdef cnp.complex128_t [:] llf_obs, u, v 1771 cdef cnp.npy_intp dim1[1] 1772 1773 u, v = zarma_innovations(endog, ar_params, ma_params, sigma2) 1774 1775 dim1[0] = nobs; 1776 llf_obs = cnp.PyArray_ZEROS(1, dim1, cnp.NPY_COMPLEX128, C) 1777 1778 const = zlog(2*NPY_PI) 1779 for i in range(nobs): 1780 llf_obs[i] = -0.5 * u[i]**2 / (sigma2 * v[i]) - 0.5 * (const + zlog(sigma2 * v[i])) 1781 1782 return np.array(llf_obs, dtype=complex) 1783