Skip to content

Commit

Permalink
REF: Remove slower versions of innovations algos.
Browse files Browse the repository at this point in the history
  • Loading branch information
@ChadFulton
ChadFulton committed Dec 27, 2018
1 parent 9aaf340 commit ab07597
Show file tree
Hide file tree
Showing 2 changed files with 0 additions and 462 deletions.
255 changes: 0 additions & 255 deletions statsmodels/tsa/innovations/_arma_innovations.pyx.in
View file
Original file line number Diff line number Diff line change
Expand Up @@ -89,109 +89,6 @@ cdef {{prefix}}toeplitz(int n, int offset0, int offset1,
out_matrix[offset0 + j, offset1 + i] = in_column[i - j]


def {{prefix}}arma_transformed_acovf({{cython_type}} [:] ar,
{{cython_type}} [:] ma,
{{cython_type}} [:] arma_acovf):
"""
arma_transformed_acovf({{cython_type}} [:] ar, {{cython_type}} [:] ma, {{cython_type}} [:] arma_acovf)

Construct the autocovariance matrix for a transformed process.

Using the autocovariance function for an ARMA process, constructs the
autocovariance matrix associated with the transformed process described
in equation (3.3.1) of _[1].

Parameters
----------
ar : ndarray
Autoregressive coefficients, including the zero lag, where the sign
convention assumes the coefficients are part of the lag polynomial on
the left-hand-side of the ARMA definition (i.e. they have the opposite
sign from the usual econometrics convention in which the coefficients
are on the right-hand-side of the ARMA definition).
ma : ndarray
Moving average coefficients, including the zero lag, where the sign
convention assumes the coefficients are part of the lag polynomial on
the right-hand-side of the ARMA definition (i.e. they have the same
sign from the usual econometrics convention in which the coefficients
are on the right-hand-side of the ARMA definition).
arma_acovf : ndarray
The vector of autocovariances of the ARMA process.

Returns
-------
acovf : ndarray
An `n` x `n` matrix containing the autocovariances of the transformed
process, where `n` is the length of the input `arma_acovf` array.

Notes
-----
The definition of this autocovariance matrix is from _[1] equation 3.3.3.

This function assumes that the variance of the ARMA innovation term is
equal to one. If this is not true, then the calling program should replace
`arma_acovf` with `arma_acovf / sigma2`, where sigma2 is that variance.

This function can be relatively slow if `arma_acovf` is very large, because
it has to allocate a matrix of that size squared. The function
`{{prefix}}arma_transformed_acovf_fast`, below, is much faster in these
cases as it only constructs a matrix corresponding to the part of the
transformed process that has time-varying autocovariances and simply
returns an array containing the autocovariance function for the remainder
of the transformed process.

References
----------
.. [1] Brockwell, Peter J., and Richard A. Davis. 2009.
Time Series: Theory and Methods. 2nd ed. 1991.
New York, NY: Springer.

"""

cdef Py_ssize_t nobs, p, q, m, m2, n, i, j, r, tmp_ix
cdef cnp.npy_intp dim1[1]
cdef cnp.npy_intp dim2[2]
cdef {{cython_type}} [:, :] acovf
cdef {{cython_type}} [:] tmp

nobs = arma_acovf.shape[0]
p = len(ar) - 1
q = len(ma) - 1
m = max(p, q)
m2 = 2 * m
n = max(m2, nobs)

dim2[0] = n;
dim2[1] = n;
acovf = cnp.PyArray_ZEROS(2, dim2, {{typenum}}, C)

# Set i, j = 1, ..., m (which is then done)
{{prefix}}toeplitz(m, 0, 0, arma_acovf, acovf)

# Set i = 1, ..., m; j = m + 1, ..., 2 m
# and i = m + 1, ..., 2 m; j = 1, ..., m
if nobs > m:
for j in range(m):
for i in range(m, m2):
acovf[i, j] = arma_acovf[i - j]
for r in range(1, p + 1):
tmp_ix = abs(r - (i - j))
acovf[i, j] = acovf[i, j] - (-ar[r] * arma_acovf[tmp_ix])
acovf[:m, m:m2] = acovf[m:m2, :m].T

# Set values for |i - j| <= q, min(i, j) = m + 1, and max(i, j) <= nobs
if nobs > m:
dim1[0] = nobs - m
tmp = cnp.PyArray_ZEROS(1, dim1, {{typenum}}, C)

for i in range(nobs - m):
for r in range(q + 1 - i):
tmp[i] = tmp[i] + ma[r] * ma[r + i]
{{prefix}}toeplitz(nobs - m, m, m, tmp, acovf)

return acovf[:nobs, :nobs]


def {{prefix}}arma_transformed_acovf_fast({{cython_type}} [:] ar,
{{cython_type}} [:] ma,
{{cython_type}} [:] arma_acovf):
Expand Down Expand Up @@ -299,86 +196,6 @@ def {{prefix}}arma_transformed_acovf_fast({{cython_type}} [:] ar,
return acovf[:n, :n], acovf2


def {{prefix}}arma_innovations_algo({{cython_type}} [:] ar_params,
{{cython_type}} [:] ma_params,
{{cython_type}} [:, :] acovf):
"""
arma_innovations_algo({{cython_type}} [:] ar_params, {{cython_type}} [:] ma_params, {{cython_type}} [:, :] acovf)

Innovations algorithm for an ARMA process.

Parameters
----------
ar_params : ndarray
Autoregressive parameters.
ma_params : ndarray
Moving average parameters.
acovf : ndarray
The autocovariance matrix of the transformed ARMA process
(see `arma_transformed_acovf`).

Returns
-------
theta : ndarray
The matrix of moving average coefficients from the innovations
algorithm.
v : ndarray
The vector of mean squared errors.

Notes
-----
The innovations algorithm is presented in _[1], section 2.5.2.

Details of the innovations algorithm applied to ARMA processes is given
in _[1] section 3.3 and in section 5.2.7.

This function can be relatively slow if `acovf` is very large. In these
cases, `arma_innovations_algo_fast` is much faster.

References
----------
.. [1] Brockwell, Peter J., and Richard A. Davis. 2009.
Time Series: Theory and Methods. 2nd ed. 1991.
New York, NY: Springer.

"""
cdef Py_ssize_t nobs, i, j, k, n, _n, m, p, q
cdef cnp.npy_intp dim1[1]
cdef cnp.npy_intp dim2[2]
cdef {{cython_type}} [:, :] theta
cdef {{cython_type}} [:] v

p = len(ar_params)
q = len(ma_params)
m = max(p, q)

nobs = acovf.shape[0] - 1

dim1[0] = nobs + 1;
v = cnp.PyArray_ZEROS(1, dim1, {{typenum}}, C)
dim2[0] = nobs + 1;
dim2[1] = nobs;
theta = cnp.PyArray_ZEROS(2, dim2, {{typenum}}, C)

v[0] = acovf[0, 0]

for n in range(nobs):
_n = n + 1
for k in range(n + 1):
if n >= m and n - k > q:
continue
theta[_n, n - k] = acovf[n + 1, k]
for j in range(k):
theta[_n, n - k] = theta[_n, n - k] - theta[k - 1 + 1, k - j - 1] * theta[_n, n - j] * v[j]
theta[_n, n - k] = theta[_n, n - k] / v[k]
v[n + 1] = acovf[n + 1, n + 1]
for i in range(n + 1):
v[n + 1] = v[n + 1] - theta[_n, n - i]**2 * v[i]

return theta, v



def {{prefix}}arma_innovations_algo_fast(int nobs,
{{cython_type}} [:] ar_params,
{{cython_type}} [:] ma_params,
Expand Down Expand Up @@ -556,78 +373,6 @@ def {{prefix}}arma_innovations_filter({{cython_type}} [:] endog,
return u


cpdef {{prefix}}arma_loglikeobs({{cython_type}} [:] endog,
{{cython_type}} [:] ar_params,
{{cython_type}} [:] ma_params,
{{cython_type}} sigma2):
"""
{{prefix}}arma_loglikeobs({{cython_type}} [:] endog, {{cython_type}} [:] ar_params, {{cython_type}} [:] ma_params, {{cython_type}} sigma2)

Calculate the loglikelihood of each observation for an ARMA process.

Parameters
----------
endog : ndarray
The observed time-series process.
ar_params : ndarray
Autoregressive parameters.
ma_params : ndarray
Moving average parameters.
sigma2 : ndarray
The ARMA innovation variance.

Returns
-------
loglike : ndarray of float
Array of loglikelihood values for each observation.

Notes
-----
Details related to computing the loglikelihood associated with an ARMA
process using the innovations algorithm are given in _[1] section 5.2.

Note that this function can be slow when the length of `endog`,
`ar_params`, or `ma_params` is large. See `arma_loglikeobs_fast` instead.

References
----------
.. [1] Brockwell, Peter J., and Richard A. Davis. 2009.
Time Series: Theory and Methods. 2nd ed. 1991.
New York, NY: Springer.

"""
cdef Py_ssize_t nobs = len(endog), i
cdef {{cython_type}} const
cdef {{cython_type}} [:] ar, ma, arma_acovf, llf_obs
cdef {{cython_type}} [:, :] acovf
cdef cnp.npy_intp dim1[1]

dim1[0] = len(ar_params) + 1;
ar = cnp.PyArray_ZEROS(1, dim1, {{typenum}}, C)
dim1[0] = len(ma_params) + 1;
ma = cnp.PyArray_ZEROS(1, dim1, {{typenum}}, C)

ar[0] = 1
for i in range(1, len(ar_params) + 1):
ar[i] = -1 * ar_params[i - 1]
ma[0] = 1
ma[1:] = ma_params

arma_acovf = arima_process.arma_acovf(ar, ma, nobs, sigma2) / sigma2
acovf = {{prefix}}arma_transformed_acovf(ar, ma, arma_acovf)
theta, v = {{prefix}}arma_innovations_algo(ar_params, ma_params, acovf)
u = {{prefix}}arma_innovations_filter(endog, ar_params, ma_params, theta)

dim1[0] = nobs;
llf_obs = cnp.PyArray_ZEROS(1, dim1, {{typenum}}, C)

const = {{combined_prefix}}log(2*NPY_PI)
for i in range(nobs):
llf_obs[i] = -0.5 * u[i]**2 / (sigma2 * v[i]) - 0.5 * (const + {{combined_prefix}}log(sigma2 * v[i]))

return np.array(llf_obs, dtype={{dtype}})


cpdef {{prefix}}arma_loglikeobs_fast({{cython_type}} [:] endog,
{{cython_type}} [:] ar_params,
{{cython_type}} [:] ma_params,
Expand Down
Loading

0 comments on commit ab07597

Please sign in to comment.