Algebra of AR, MA models
Contents
Algebra of AR, MA models¶
import warnings
warnings.simplefilter("ignore", FutureWarning)
import matplotlib.pyplot as plt
import numpy as np
from statsmodels.tsa.arima_process import arma_acovf, ArmaProcess
from scipy.linalg import toeplitz
from matplotlib import rc
rc('font',**{'family':'serif','serif':['Palatino']})
rc('text', usetex=True)
plt.rcParams.update({'ytick.left' : True,
"xtick.bottom" : True,
"ytick.major.size": 0,
"ytick.major.width": 0,
"xtick.major.size": 0,
"xtick.major.width": 0,
"ytick.direction": 'in',
"xtick.direction": 'in',
'ytick.major.right': False,
'xtick.major.top': True,
'xtick.top': True,
'ytick.right': True,
'ytick.labelsize': 18,
'xtick.labelsize': 18
})
np.random.seed(12345)
MA models¶
Example:
\[z_t = \varepsilon_t + 0.4 \varepsilon_{t-1}\]
define in Python
theta = np.array([.4])
ar = np.r_[1] # coefficient on z_t
ma = np.r_[1, theta] # coefficients on e(t) and e(t-1)
ma1_process = ArmaProcess(ar, ma)
print(ma1_process)
ArmaProcess
AR: [1.0]
MA: [1.0, 0.4]
# what is ma1_process?
type(ma1_process)
statsmodels.tsa.arima_process.ArmaProcess
# what are ma1_process's attributes?
[x for x in dir(ma1_process) if not x.startswith('__')]
['acf',
'acovf',
'ar',
'arcoefs',
'arma2ar',
'arma2ma',
'arpoly',
'arroots',
'from_coeffs',
'from_estimation',
'from_roots',
'generate_sample',
'impulse_response',
'invertroots',
'isinvertible',
'isstationary',
'ma',
'macoefs',
'mapoly',
'maroots',
'nobs',
'pacf',
'periodogram']
ma1_process.mapoly # MA polynomial
\[x \mapsto \text{1.0} + \text{0.4}\,x\]
ma1_process.maroots # MA roots
array([-2.5])
the MA root(s) solve the MA poly(nomial)
ma_root = -2.5
1.0 + 0.4 * ma_root
0.0
Consider a MA(2) model
theta = np.array([.6, -.3]) # coefficients on e(t-1) and e(t-2)
ar = np.r_[1] # coefficient on z_t
ma = np.r_[1, theta]
ma2_process = ArmaProcess(ar, ma)
ma2_process.mapoly
\[x \mapsto \text{1.0} + \text{0.6}\,x - \text{0.3}\,x^{2}\]
ma2_process.maroots
array([-1.081666, 3.081666])
ma_root1 = ma2_process.maroots[0]
expression1 = 1.0 + (0.6) * (ma_root1) - (0.3)*(ma_root1)**2
expression1
1.6653345369377348e-16
An aside:
Numerically, it is often difficult to distinguish very small numbers from exectly 0. We can confirm that the value is almost zero using the following assertion:
# using numpy
from numpy.testing import assert_almost_equal
assert_almost_equal(0,expression1, err_msg="expression1 is NOT almost equal to 0!")
# here is what happens when the assertion is false
assert_almost_equal(0.01,expression1, err_msg="expression1 is NOT almost equal to 0!")
---------------------------------------------------------------------------
AssertionError Traceback (most recent call last)
Input In [69], in <module>
1 # here is what happens when the assertion is false
----> 2 assert_almost_equal(0.01,expression1, err_msg="expression1 is NOT almost equal to 0!")
File ~\AppData\Local\Continuum\anaconda3\envs\teaching\lib\site-packages\numpy\testing\_private\utils.py:599, in assert_almost_equal(actual, desired, decimal, err_msg, verbose)
597 pass
598 if abs(desired - actual) >= 1.5 * 10.0**(-decimal):
--> 599 raise AssertionError(_build_err_msg())
AssertionError:
Arrays are not almost equal to 7 decimals expression1 is NOT almost equal to 0!
ACTUAL: 0.01
DESIRED: 4.440892098500626e-16
# using standard python
very_small_number = 0.00000001
assert expression1 < very_small_number, "expression1 is NOT larger than the 'very_small_number'!"
# here is what happens when the assertion is false
very_small_number = 0.00000001
assert expression1 > very_small_number, "expression1 is NOT larger than the 'very_small_number'!"
---------------------------------------------------------------------------
AssertionError Traceback (most recent call last)
Input In [71], in <module>
1 # here is what happens when the assertion is false
2 very_small_number = 0.00000001
----> 3 assert expression1 > very_small_number, "expression1 is NOT larger than the 'very_small_number'!"
AssertionError: expression1 is NOT larger than the 'very_small_number'!
ma_root2 = ma2_process.maroots[1]
expression2 = 1.0 + (0.6) * (ma_root2) - (0.3)*(ma_root2)**2
expression2
0.0
AR models¶
Example:
\[z_t = 0.8 z_{t-1} + \varepsilon_t\]
Define in Python
alpha = .8 # coefficient on z(t-1)
ar = np.r_[1, -alpha]
ma = np.r_[1] # coefficient on e(t)
ar1_process = ArmaProcess(ar, ma)
print(ar1_process)
ArmaProcess
AR: [1.0, -0.8]
MA: [1.0]
ar1_process.arpoly
\[x \mapsto \text{1.0} - \text{0.8}\,x\]
ar1_process.arroots
array([1.25])
the AR root(s) solves the AR poly(nomial)
ar_root = 1.25
1.0 - 0.8 * (ar_root)
0.0
alpha = np.array([.4, .1])
ar = np.r_[1, -alpha]
ma = np.r_[1]
ar2_process = ArmaProcess(ar, ma)
ar2_process.arpoly
\[x \mapsto \text{1.0} - \text{0.4}\,x - \text{0.1}\,x^{2}\]
ar2_process.arroots
array([-5.74165739, 1.74165739])
ar_root1 = ar2_process.arroots[0]
expression1 = 1.0 - (0.4) * (ar_root1) - (0.1)*(ar_root1)**2
expression1
4.440892098500626e-16
ar_root2 = ar2_process.arroots[1]
expression2 = 1.0 - (0.4) * (ar_root2) - (0.1)*(ar_root2)**2
expression2
-1.6653345369377348e-16
Lag operator¶
\[ L z_{t} = z_{t-1} \]
\[ L^2 z_{t} = L(L z_t) = Lz_{t-1} = z_{t-2} \]
\[ L^n z_{t} = z_{t-n} \]
Note
aka Backshift operator
polynomial of order \(p\) in the lag operator
\[ \alpha(L) = \alpha_0 + \alpha_1 L + \alpha_2 L^2 + \cdots + \alpha_p L^p\]
weighted moving average of \(z_t\)
\[\begin{split}
\begin{align}
\alpha(L)z_t &= z_t(\alpha_0 + \alpha_1 L + \alpha_2 L^2 + \cdots + \alpha_p L^p) \\
&= \alpha_0 z_t + \alpha_1 z_{t-1} + \alpha_2 z_{t-2} + \cdots + \alpha_p z_{t-p}
\end{align}
\end{split}\]
AR models in lag operator notation¶
AR(1)¶
\[
z_{t} = \alpha z_{t-1} + \varepsilon_{t}
\]
\[
(1- \alpha L) z_{t} = \varepsilon_{t}
\]
\[
\text{if } |\alpha|< 1 \rightarrow z_{t} = \frac{1}{(1- \alpha L)}\varepsilon_{t}
\]
geometric series result
\[ \sum_{j=0}^{\infty} x^j = 1 + x + x^2+\cdots= \frac{1}{1-x}, \;\;\; \text{if } |x|<1\]
in our case \(x = \alpha L\)
\[ \frac{1}{(1- \alpha L)} = \sum_{j=0}^{\infty} (\alpha L)^j = 1 + \alpha L + \alpha^2 L^2+\cdots\]
\[
z_{t} = \frac{1}{(1- \alpha L)}\varepsilon_{t} = \varepsilon_{t} + \alpha \varepsilon_{t-1} + \alpha^2 \varepsilon_{t-2}+\cdots = \sum_{i=0}^{\infty} \alpha^i \varepsilon_{t-i}
\]
AR(1) is a MA(\(\infty\))
another derivation of the variance of AR(1)
\[
z_{t} = \varepsilon_{t} + \alpha \varepsilon_{t-1} + \alpha^2 \varepsilon_{t-2}+\cdots
\]
\[
\operatorname{var}(z_{t}) = \operatorname{var}(\varepsilon_{t}) + \alpha^2 \operatorname{var}(\varepsilon_{t-1}) + \alpha^4 \operatorname{var}(\varepsilon_{t-2}) + \cdots = \sigma^2 \sum_{i=0}^{\infty} (\alpha^2)^i = \frac{\sigma^2}{1-\alpha^2}
\]
AR(1) is stationary if \(|\alpha|<1\)
equivalent to polynomilal
\[\alpha(x) = 1 - \alpha x\]
having all roots \(|x|>1\) (outside the unit circle)
ar1_process.arpoly
\[x \mapsto \text{1.0} - \text{0.8}\,x\]
ar1_process.arroots
array([1.25])
np.abs(ar1_process.arroots)>1
array([ True])
print(ar1_process.isstationary)
True
AR(q) process is stationary if all roots \(z\) of the polynomial
\[ \alpha(x) = 1 - \alpha_1 x - \alpha_2 x^2 - \cdots - \alpha_q x^q\]
are outside the unit circle (\(|x|>1\))
ar2_process.arpoly
\[x \mapsto \text{1.0} - \text{0.4}\,x - \text{0.1}\,x^{2}\]
ar2_process.arroots
array([-5.74165739, 1.74165739])
all(np.abs(ar2_process.arroots)>1)
True
print(ar2_process.isstationary)
True
MA models in lag operator notation¶
MA(1)¶
\[
z_{t} = \varepsilon_{t} + \theta \varepsilon_{t-1}
\]
\[
z_{t} = (1 + \theta L) \varepsilon_{t}
\]
\[
\varepsilon_{t} = \frac{1}{(1 + \theta L)} z_{t} = z_{t} - \theta z_{t-1} + \theta^2 z_{t-2}+\cdots = \sum_{i=0}^{\infty} (-\theta)^i z_{t-i}, \;\;\; \textbf{if } |\theta|<1
\]
if \(|\theta|<1\) then \(\varepsilon_t\) can be represented in terms of the current and past values of \(z_t\)
the process is invertible
MA(1) process for \(z_t\) is a stable AR(1) process for \(\varepsilon\)
MA models are unidentified¶
MA(1)
model 1
\[\begin{split}
\begin{align}
z_{t} = \varepsilon_{t} &+ \theta \varepsilon_{t-1} , \; \operatorname{var} (\varepsilon_{t})=\sigma^2 \\
&\operatorname{var} (z) = \sigma^2 + \theta^2 \sigma^2\\
&\operatorname{cov}(z_t, z_{t-1}) = \theta \sigma^2\\
&\operatorname{cov}(z_t, z_{t-h}) = 0
\end{align}
\end{split}\]
model 2
\[
\hat{z}_{t} = \hat{\varepsilon}_{t} + \frac{1}{\theta} \hat{\varepsilon}_{t-1}, \;\;\; \operatorname{var} (\hat{\varepsilon}_{t})=\sigma^2 \theta^2
\]
equivalent to:
\[
\hat{z}_{t} = \theta \varepsilon_{t} + \varepsilon_{t-1}, \;\;\; \operatorname{var} (\varepsilon_{t})=\sigma^2
\]
\[
\operatorname{var} (\hat{z}) = \theta^2 \sigma^2 + \frac{1}{\theta^2}\theta^2\sigma^2 = \theta^2 \sigma^2 + \sigma^2 =\operatorname{var} (z)
\]
\[
\operatorname{cov}(\hat{z}_t, \hat{z}_{t-1}) = \frac{1}{\theta} \sigma^2 \theta^2 = \sigma^2 \theta = \operatorname{cov}(z_t, z_{t-1})
\]
model 1 and model 2 are indistinguishable from variances and autocovariances
two different values of the parameters (\(\theta\) and \(\sigma\)) produce the same value of the likelihood
if model 1 is invertable (\(|\theta|<1\)), model 2 is not (\(\frac{1}{|\theta|}>1\))
the invertible condition may be desirable to impose: think of \(\varepsilon_t\) as an economic shock (e.g. TFP). From the data (\(z\)) we observe only the effect of the shock, not the shock itself. Inverability implies that the values of the shock can be recovered from past values of \(z_t\)
ma1_process
ArmaProcess([1.0], [1.0, 0.4], nobs=100) at 0x2297a713460
ma1_process.isinvertible
True
Here is an example of a non-invertible MA(1) process
non_invertible_ma1_process = ArmaProcess(1, [0.4, 1], )
non_invertible_ma1_process
ArmaProcess([1.0], [0.4, 1.0], nobs=100) at 0x2297e62ff10
Compare the covariances with the one of the MA(1) process defined above:
non_invertible_ma1_process.acovf(3)
array([1.16, 0.4 , 0. ])
ma1_process.acovf(3)
array([1.16, 0.4 , 0. ])
The non-invertible process is not invertible
non_invertible_ma1_process.isinvertible
False
For a AR(2) models with coefficients \(\alpha_1\) and \(\alpha_2\), the stability conditions are the following: $$
(2)¶\[\begin{align}
\alpha_1 + \alpha_2 &< 1\\
\alpha_2 - \alpha_1 &< 1\\
-1<\alpha_2 &<1
\end{align}\]
$$
This is known as the stability triangle.
# stability triangle
ar2_process.arcoefs.sum()<1, ar2_process.arcoefs[1] - ar2_process.arcoefs[0]<1, np.abs(ar2_process.arcoefs[1])<1
(True, True, True)