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{{py: 12 13TYPES = { 14 "s": ("cnp.float32_t", "np.float32", "cnp.NPY_FLOAT32"), 15 "d": ("cnp.float64_t", "float", "cnp.NPY_FLOAT64"), 16 "c": ("cnp.complex64_t", "np.complex64", "cnp.NPY_COMPLEX64"), 17 "z": ("cnp.complex128_t", "complex", "cnp.NPY_COMPLEX128"), 18} 19 20}} 21 22# Typical imports 23import numpy as np 24import warnings 25from statsmodels.tsa import arima_process 26from statsmodels.src.math cimport NPY_PI, dlog, zlog 27cimport numpy as cnp 28cimport cython 29 30cnp.import_array() 31 32cdef int C = 0 33 34 35{{for prefix, types in TYPES.items()}} 36{{py:cython_type, dtype, typenum = types}} 37{{py: 38combined_prefix = prefix 39combined_cython_type = cython_type 40if prefix == 'c': 41 combined_prefix = 'z' 42 combined_cython_type = 'np.complex128_t' 43if prefix == 's': 44 combined_prefix = 'd' 45 combined_cython_type = 'np.float64_t' 46}} 47 48 49cdef {{prefix}}toeplitz(int n, int offset0, int offset1, 50 {{cython_type}} [:] in_column, 51 {{cython_type}} [:, :] out_matrix): 52 """ 53 toeplitz(int n, int o0, int o1, {{cython_type}} [:] in_column, {{cython_type}} [:, :] out_matrix 54 55 Construct a Toeplitz block in a matrix in place 56 57 Parameters 58 ---------- 59 n : int 60 The number of entries of `in_column` to use. 61 offset0 : int 62 The row offset for `out_matrix` at which to begin writing the Toeplitz 63 block. 64 offset1 : int 65 The column offset for `out_matrix` at which to begin writing the 66 Toeplitz block. 67 in_column : ndarray 68 The column used to construct the Toeplitz block. 69 out_matrix : ndarray 70 The matrix in which to write a Toeplitz block 71 72 Notes 73 ----- 74 This function constructs the Toeplitz block in-place and does not return 75 any output. 76 """ 77 cdef Py_ssize_t i, j 78 79 for i in range(n): 80 for j in range(i + 1): 81 out_matrix[offset0 + i, offset1 + j] = in_column[i - j] 82 if i != j: 83 # Note: scipy by default does complex conjugate, but not 84 # necessary here since we're only dealing with covariances, 85 # which will be real (except in the case of complex-step 86 # differentiation, but in that case we do not want to apply 87 # the conjugation anyway) 88 out_matrix[offset0 + j, offset1 + i] = in_column[i - j] 89 90 91cpdef {{prefix}}arma_transformed_acovf_fast({{cython_type}} [:] ar, 92 {{cython_type}} [:] ma, 93 {{cython_type}} [:] arma_acovf): 94 """ 95 arma_transformed_acovf_fast({{cython_type}} [:] ar, {{cython_type}} [:] ma, {{cython_type}} [:] arma_acovf) 96 97 Quickly construct the autocovariance matrix for a transformed process. 98 99 Using the autocovariance function for an ARMA process, constructs the 100 autocovariances associated with the transformed process described 101 in equation (3.3.1) of _[1] in a memory efficient, and so fast, way. 102 103 Parameters 104 ---------- 105 ar : ndarray 106 Autoregressive coefficients, including the zero lag, where the sign 107 convention assumes the coefficients are part of the lag polynomial on 108 the left-hand-side of the ARMA definition (i.e. they have the opposite 109 sign from the usual econometrics convention in which the coefficients 110 are on the right-hand-side of the ARMA definition). 111 ma : ndarray 112 Moving average coefficients, including the zero lag, where the sign 113 convention assumes the coefficients are part of the lag polynomial on 114 the right-hand-side of the ARMA definition (i.e. they have the same 115 sign from the usual econometrics convention in which the coefficients 116 are on the right-hand-side of the ARMA definition). 117 arma_acovf : ndarray 118 The vector of autocovariances of the ARMA process. 119 120 Returns 121 ------- 122 acovf : ndarray 123 A matrix containing the autocovariances of the portion of the 124 transformed process with time-varying autocovariances. Its dimension 125 is equal to `min(m * 2, n)` where `m = max(len(ar) - 1, len(ma) - 1)` 126 and `n` is the length of the input array `arma_acovf`. It is important 127 to note that only the values in the first `m` columns or `m` rows are 128 valid. In particular, the entries in the block `acovf[m:, m:]` should 129 not be used in any case (and in fact will always be equal to zeros). 130 acovf2 : ndarray 131 An array containing the autocovariance function of the portion of the 132 transformed process with time-invariant autocovariances. Its dimension 133 is equal to `max(n - m, 0)` where `n` is the length of the input 134 array `arma_acovf`. 135 136 Notes 137 ----- 138 The definition of this autocovariance matrix is from _[1] equation 3.3.3. 139 140 This function assumes that the variance of the ARMA innovation term is 141 equal to one. If this is not true, then the calling program should replace 142 `arma_acovf` with `arma_acovf / sigma2`, where sigma2 is that variance. 143 144 This function is relatively fast even when `arma_acovf` is large, since 145 it only constructs the full autocovariance matrix for a generally small 146 subset of observations. The trade-off is that the output of this function 147 is somewhat more difficult to use. 148 149 References 150 ---------- 151 .. [1] Brockwell, Peter J., and Richard A. Davis. 2009. 152 Time Series: Theory and Methods. 2nd ed. 1991. 153 New York, NY: Springer. 154 """ 155 cdef Py_ssize_t nobs, p, q, m, m2, n, i, j, r, tmp_ix 156 cdef cnp.npy_intp dim1[1] 157 cdef cnp.npy_intp dim2[2] 158 cdef {{cython_type}} [:, :] acovf 159 cdef {{cython_type}} [:] acovf2 160 161 nobs = arma_acovf.shape[0] 162 p = len(ar) - 1 163 q = len(ma) - 1 164 m = max(p, q) 165 m2 = 2 * m 166 n = min(m2, nobs) 167 168 dim2[0] = m2; 169 dim2[1] = m2; 170 acovf = cnp.PyArray_ZEROS(2, dim2, {{typenum}}, C) 171 dim1[0] = max(nobs - m, 0); 172 acovf2 = cnp.PyArray_ZEROS(1, dim1, {{typenum}}, C) 173 174 # Set i, j = 1, ..., m (which is then done) 175 {{prefix}}toeplitz(m, 0, 0, arma_acovf, acovf) 176 177 # Set i = 1, ..., m; j = m + 1, ..., 2 m 178 # and i = m + 1, ..., 2 m; j = 1, ..., m 179 if nobs > m: 180 for j in range(m): 181 for i in range(m, m2): 182 acovf[i, j] = arma_acovf[i - j] 183 for r in range(1, p + 1): 184 tmp_ix = abs(r - (i - j)) 185 acovf[i, j] = acovf[i, j] - (-ar[r] * arma_acovf[tmp_ix]) 186 acovf[:m, m:m2] = acovf[m:m2, :m].T 187 188 # Set values for |i - j| <= q, min(i, j) = m + 1, and max(i, j) <= nobs 189 if nobs > m: 190 for i in range(nobs - m): 191 for r in range(q + 1 - i): 192 acovf2[i] = acovf2[i] + ma[r] * ma[r + i] 193 194 return acovf[:n, :n], acovf2 195 196 197cpdef {{prefix}}arma_innovations_algo_fast(int nobs, 198 {{cython_type}} [:] ar_params, 199 {{cython_type}} [:] ma_params, 200 {{cython_type}} [:, :] acovf, 201 {{cython_type}} [:] acovf2): 202 """ 203 arma_innovations_algo_fast(int nobs, {{cython_type}} [:] ar_params, {{cython_type}} [:] ma_params, {{cython_type}} [:, :] acovf, {{cython_type}} [:] acovf2) 204 205 Quickly apply innovations algorithm for an ARMA process. 206 207 Parameters 208 ---------- 209 ar_params : ndarray 210 Autoregressive parameters. 211 ma_params : ndarray 212 Moving average parameters. 213 acovf : ndarray 214 An `m * 2` x `m * 2` autocovariance matrix at least the first `m` 215 columns filled in, where `m = max(len(ar_params), ma_params)` 216 (see `arma_transformed_acovf_fast`). 217 acovf2 : ndarray 218 A `max(0, nobs - m)` length vector containing the autocovariance 219 function associated with the final `nobs - m` observations. 220 221 Returns 222 ------- 223 theta : ndarray 224 The matrix of moving average coefficients from the innovations 225 algorithm. 226 v : ndarray 227 The vector of mean squared errors. 228 229 Notes 230 ----- 231 The innovations algorithm is presented in _[1], section 2.5.2. 232 233 Details of the innovations algorithm applied to ARMA processes is given 234 in _[1] section 3.3 and in section 5.2.7. 235 236 This function is relatively fast even with a large number of observations 237 since we can exploit a number of known zeros in the theta array. 238 239 References 240 ---------- 241 .. [1] Brockwell, Peter J., and Richard A. Davis. 2009. 242 Time Series: Theory and Methods. 2nd ed. 1991. 243 New York, NY: Springer. 244 """ 245 cdef Py_ssize_t i, j, k, n, _n, m, m2, p, q, start, start2 246 cdef cnp.npy_intp dim1[1] 247 cdef cnp.npy_intp dim2[2] 248 cdef {{cython_type}} [:, :] theta 249 cdef {{cython_type}} [:] v 250 251 p = len(ar_params) 252 q = len(ma_params) 253 m = max(p, q) 254 m2 = 2 * m 255 256 dim1[0] = nobs; 257 v = cnp.PyArray_ZEROS(1, dim1, {{typenum}}, C) 258 dim2[0] = nobs; 259 dim2[1] = m + 1; 260 theta = cnp.PyArray_ZEROS(2, dim2, {{typenum}}, C) 261 262 if m > 0: 263 v[0] = acovf[0, 0] 264 else: 265 # (handle the corner case where p = q = 0) 266 v[0] = acovf2[0] 267 268 for n in range(nobs - 1): 269 _n = n + 1 270 271 start = 0 if n < m else n + 1 - q 272 for k in range(start, n + 1): 273 # See Brockwell and Davis, p. 100-101 274 # (here we have weak rather than strict inequality due to Python's 275 # zero indexing) 276 if n >= m and n - k >= q: 277 continue 278 279 if n + 1 < m2 and k < m: 280 theta[_n, n - k] = acovf[n + 1, k] 281 else: 282 theta[_n, n - k] = acovf2[n + 1 - k] 283 284 start2 = 0 if n < m else n - m 285 for j in range(start2, k): 286 if n - j < m + 1: 287 theta[_n, n - k] = theta[_n, n - k] - theta[k - 1 + 1, k - j - 1] * theta[_n, n - j] * v[j] 288 theta[_n, n - k] = theta[_n, n - k] / v[k] 289 290 if n + 1 < m: 291 v[n + 1] = acovf[n + 1, n + 1] 292 else: 293 v[n + 1] = acovf2[0] 294 start = max(0, n - (m + 1) + 2) 295 for i in range(start, n + 1): 296 v[n + 1] = v[n + 1] - theta[_n, n - i]**2 * v[i] 297 298 return theta, v 299 300 301cpdef {{prefix}}arma_innovations_filter({{cython_type}} [:] endog, 302 {{cython_type}} [:] ar_params, 303 {{cython_type}} [:] ma_params, 304 {{cython_type}} [:, :] theta): 305 """ 306 arma_innovations_filter({{cython_type}} [:] endog, {{cython_type}} [:] ar_params, {{cython_type}} [:] ma_params, {{cython_type}} [:, :] theta): 307 308 Innovations filter for an ARMA process. 309 310 Parameters 311 ---------- 312 endog : ndarray 313 The observed time-series process. 314 ar_params : ndarray 315 Autoregressive parameters. 316 ma_params : ndarray 317 Moving average parameters. 318 theta : ndarray 319 The matrix of moving average coefficients from the innovations 320 algorithm (see `arma_innovations_algo` or `arma_innovations_algo_fast`) 321 322 Returns 323 ------- 324 u : ndarray 325 The vector of innovations: the one-step-ahead prediction errors from 326 applying the innovations algorithm. 327 328 Notes 329 ----- 330 The innovations algorithm is presented in _[1], section 2.5.2. 331 332 Details of the innovations algorithm applied to ARMA processes is given 333 in _[1] section 3.3 and in section 5.2.7. 334 335 References 336 ---------- 337 .. [1] Brockwell, Peter J., and Richard A. Davis. 2009. 338 Time Series: Theory and Methods. 2nd ed. 1991. 339 New York, NY: Springer. 340 """ 341 cdef Py_ssize_t nobs, i, k, j, m, p, q 342 cdef cnp.npy_intp dim1[1] 343 cdef {{cython_type}} [:] u 344 cdef {{cython_type}} hat 345 346 p = len(ar_params) 347 q = len(ma_params) 348 m = max(p, q) 349 350 nobs = theta.shape[0] 351 k = theta.shape[1] 352 353 dim1[0] = nobs; 354 u = cnp.PyArray_ZEROS(1, dim1, {{typenum}}, C) 355 u[0] = endog[0] 356 357 for i in range(1, nobs): 358 hat = 0 359 if i < m: 360 for j in range(i): 361 hat = hat + theta[i, j] * u[i - j - 1] 362 else: 363 for j in range(p): 364 hat = hat + ar_params[j] * endog[i - j - 1] 365 for j in range(q): 366 hat = hat + theta[i, j] * u[i - j - 1] 367 u[i] = endog[i] - hat 368 369 return u 370 371 372cpdef {{prefix}}arma_innovations({{cython_type}} [:] endog, 373 {{cython_type}} [:] ar_params, 374 {{cython_type}} [:] ma_params, 375 {{cython_type}} sigma2): 376 """ 377 arma_innovations({{cython_type}} [:] endog, {{cython_type}} [:] ar_params, {{cython_type}} [:] ma_params): 378 379 Compute innovations and variances based on an ARMA process. 380 381 Parameters 382 ---------- 383 endog : ndarray 384 The observed time-series process. 385 ar_params : ndarray 386 Autoregressive parameters. 387 ma_params : ndarray 388 Moving average parameters. 389 sigma2 : ndarray 390 The ARMA innovation variance. 391 392 Returns 393 ------- 394 u : ndarray 395 The vector of innovations: the one-step-ahead prediction errors from 396 applying the innovations algorithm. 397 v : ndarray 398 The vector of innovation variances. 399 400 Notes 401 ----- 402 The innovations algorithm is presented in _[1], section 2.5.2. 403 404 References 405 ---------- 406 .. [1] Brockwell, Peter J., and Richard A. Davis. 2009. 407 Time Series: Theory and Methods. 2nd ed. 1991. 408 New York, NY: Springer. 409 """ 410 cdef Py_ssize_t nobs = len(endog), i 411 cdef {{cython_type}} const 412 cdef {{cython_type}} [:] ar, ma, arma_acovf, llf_obs, acovf2, u, v 413 cdef {{cython_type}} [:, :] acovf 414 cdef cnp.npy_intp dim1[1] 415 416 dim1[0] = len(ar_params) + 1; 417 ar = cnp.PyArray_ZEROS(1, dim1, {{typenum}}, C) 418 dim1[0] = len(ma_params) + 1; 419 ma = cnp.PyArray_ZEROS(1, dim1, {{typenum}}, C) 420 421 ar[0] = 1 422 for i in range(1, len(ar_params) + 1): 423 ar[i] = -1 * ar_params[i - 1] 424 ma[0] = 1 425 ma[1:] = ma_params 426 427 arma_acovf = arima_process.arma_acovf(ar, ma, nobs, sigma2, dtype={{dtype}}) / sigma2 428 acovf, acovf2 = {{prefix}}arma_transformed_acovf_fast(ar, ma, arma_acovf) 429 theta, v = {{prefix}}arma_innovations_algo_fast(nobs, ar_params, ma_params, acovf, acovf2) 430 u = {{prefix}}arma_innovations_filter(endog, ar_params, ma_params, theta) 431 432 return u, v 433 434 435cpdef {{prefix}}arma_loglikeobs_fast({{cython_type}} [:] endog, 436 {{cython_type}} [:] ar_params, 437 {{cython_type}} [:] ma_params, 438 {{cython_type}} sigma2): 439 """ 440 {{prefix}}arma_loglikeobs_fast({{cython_type}} [:] endog, {{cython_type}} [:] ar_params, {{cython_type}} [:] ma_params, {{cython_type}} sigma2) 441 442 Quickly calculate the loglikelihood of each observation for an ARMA process 443 444 Parameters 445 ---------- 446 endog : ndarray 447 The observed time-series process. 448 ar_params : ndarray 449 Autoregressive parameters. 450 ma_params : ndarray 451 Moving average parameters. 452 sigma2 : ndarray 453 The ARMA innovation variance. 454 455 Returns 456 ------- 457 loglike : ndarray of float 458 Array of loglikelihood values for each observation. 459 460 Notes 461 ----- 462 Details related to computing the loglikelihood associated with an ARMA 463 process using the innovations algorithm are given in _[1] section 5.2. 464 465 References 466 ---------- 467 .. [1] Brockwell, Peter J., and Richard A. Davis. 2009. 468 Time Series: Theory and Methods. 2nd ed. 1991. 469 New York, NY: Springer. 470 """ 471 472 cdef Py_ssize_t nobs = len(endog), i 473 cdef {{cython_type}} const 474 cdef {{cython_type}} [:] llf_obs, u, v 475 cdef cnp.npy_intp dim1[1] 476 477 u, v = {{prefix}}arma_innovations(endog, ar_params, ma_params, sigma2) 478 479 dim1[0] = nobs; 480 llf_obs = cnp.PyArray_ZEROS(1, dim1, {{typenum}}, C) 481 482 const = {{combined_prefix}}log(2*NPY_PI) 483 for i in range(nobs): 484 llf_obs[i] = -0.5 * u[i]**2 / (sigma2 * v[i]) - 0.5 * (const + {{combined_prefix}}log(sigma2 * v[i])) 485 486 return np.array(llf_obs, dtype={{dtype}}) 487 488{{endfor}} 489