Source code for quantecon.arma

Provides functions for working with and visualizing scalar ARMA processes.

TODO: 1. Fix warnings concerning casting complex variables back to floats

import numpy as np
from numpy import conj, pi
from .util import check_random_state

[docs]class ARMA: r""" This class represents scalar ARMA(p, q) processes. If phi and theta are scalars, then the model is understood to be .. math:: X_t = \phi X_{t-1} + \epsilon_t + \theta \epsilon_{t-1} where :math:`\epsilon_t` is a white noise process with standard deviation :math:`\sigma`. If phi and theta are arrays or sequences, then the interpretation is the ARMA(p, q) model .. math:: X_t = \phi_1 X_{t-1} + ... + \phi_p X_{t-p} + \epsilon_t + \theta_1 \epsilon_{t-1} + ... + \theta_q \epsilon_{t-q} where * :math:`\phi = (\phi_1, \phi_2,..., \phi_p)` * :math:`\theta = (\theta_1, \theta_2,..., \theta_q)` * :math:`\sigma` is a scalar, the standard deviation of the white noise Parameters ---------- phi : scalar or iterable or array_like(float) Autocorrelation values for the autocorrelated variable. See above for explanation. theta : scalar or iterable or array_like(float) Autocorrelation values for the white noise of the model. See above for explanation sigma : scalar(float) The standard deviation of the white noise Attributes ---------- phi, theta, sigma : see Parmeters ar_poly : array_like(float) The polynomial form that is needed by scipy.signal to do the processing we desire. Corresponds with the phi values ma_poly : array_like(float) The polynomial form that is needed by scipy.signal to do the processing we desire. Corresponds with the theta values """ def __init__(self, phi, theta=0, sigma=1): self._phi, self._theta = phi, theta self.sigma = sigma self.set_params() def __repr__(self): m = "ARMA(phi=%s, theta=%s, sigma=%s)" return m % (self.phi, self.theta, self.sigma) def __str__(self): m = "An ARMA({p}, {q}) process" p = np.asarray(self.phi).size q = np.asarray(self.theta).size return m.format(p=p, q=q) # Special latex print method for working in notebook def _repr_latex_(self): m = r"$X_t = " phi = np.atleast_1d(self.phi) theta = np.atleast_1d(self.theta) rhs = "" for (tm, phi_p) in enumerate(phi): # don't include terms if they are equal to zero if abs(phi_p) > 1e-12: rhs += r"%+g X_{t-%i}" % (phi_p, tm+1) if rhs[0] == "+": rhs = rhs[1:] # remove initial `+` if phi_1 was positive rhs += r" + \epsilon_t" for (tm, th_q) in enumerate(theta): # don't include terms if they are equal to zero if abs(th_q) > 1e-12: rhs += r"%+g \epsilon_{t-%i}" % (th_q, tm+1) return m + rhs + "$" @property def phi(self): return self._phi @phi.setter def phi(self, new_value): self._phi = new_value self.set_params() @property def theta(self): return self._theta @theta.setter def theta(self, new_value): self._theta = new_value self.set_params()
[docs] def set_params(self): r""" Internally, scipy.signal works with systems of the form .. math:: ar_{poly}(L) X_t = ma_{poly}(L) \epsilon_t where L is the lag operator. To match this, we set .. math:: ar_{poly} = (1, -\phi_1, -\phi_2,..., -\phi_p) ma_{poly} = (1, \theta_1, \theta_2,..., \theta_q) In addition, ar_poly must be at least as long as ma_poly. This can be achieved by padding it out with zeros when required. """ # === set up ma_poly === # ma_poly = np.asarray(self._theta) self.ma_poly = np.insert(ma_poly, 0, 1) # The array (1, theta) # === set up ar_poly === # if np.isscalar(self._phi): ar_poly = np.array(-self._phi) else: ar_poly = -np.asarray(self._phi) self.ar_poly = np.insert(ar_poly, 0, 1) # The array (1, -phi) # === pad ar_poly with zeros if required === # if len(self.ar_poly) < len(self.ma_poly): temp = np.zeros(len(self.ma_poly) - len(self.ar_poly)) self.ar_poly = np.hstack((self.ar_poly, temp))
[docs] def impulse_response(self, impulse_length=30): """ Get the impulse response corresponding to our model. Returns ------- psi : array_like(float) psi[j] is the response at lag j of the impulse response. We take psi[0] as unity. """ from scipy.signal import dimpulse sys = self.ma_poly, self.ar_poly, 1 times, psi = dimpulse(sys, n=impulse_length) psi = psi[0].flatten() # Simplify return value into flat array return psi
[docs] def spectral_density(self, two_pi=True, res=1200): r""" Compute the spectral density function. The spectral density is the discrete time Fourier transform of the autocovariance function. In particular, .. math:: f(w) = \sum_k \gamma(k) \exp(-ikw) where gamma is the autocovariance function and the sum is over the set of all integers. Parameters ---------- two_pi : Boolean, optional Compute the spectral density function over :math:`[0, \pi]` if two_pi is False and :math:`[0, 2 \pi]` otherwise. Default value is True res : scalar or array_like(int), optional(default=1200) If res is a scalar then the spectral density is computed at `res` frequencies evenly spaced around the unit circle, but if res is an array then the function computes the response at the frequencies given by the array Returns ------- w : array_like(float) The normalized frequencies at which h was computed, in radians/sample spect : array_like(float) The frequency response """ from scipy.signal import freqz w, h = freqz(self.ma_poly, self.ar_poly, worN=res, whole=two_pi) spect = h * conj(h) * self.sigma**2 return w, spect
[docs] def autocovariance(self, num_autocov=16): """ Compute the autocovariance function from the ARMA parameters over the integers range(num_autocov) using the spectral density and the inverse Fourier transform. Parameters ---------- num_autocov : scalar(int), optional(default=16) The number of autocovariances to calculate """ spect = self.spectral_density()[1] acov = np.fft.ifft(spect).real # num_autocov should be <= len(acov) / 2 return acov[:num_autocov]
[docs] def simulation(self, ts_length=90, random_state=None): """ Compute a simulated sample path assuming Gaussian shocks. Parameters ---------- ts_length : scalar(int), optional(default=90) Number of periods to simulate for random_state : int or np.random.RandomState, optional Random seed (integer) or np.random.RandomState instance to set the initial state of the random number generator for reproducibility. If None, a randomly initialized RandomState is used. Returns ------- vals : array_like(float) A simulation of the model that corresponds to this class """ from scipy.signal import dlsim random_state = check_random_state(random_state) sys = self.ma_poly, self.ar_poly, 1 u = random_state.randn(ts_length, 1) * self.sigma vals = dlsim(sys, u)[1] return vals.flatten()