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