Objective
Impulse response functions (IRFs) are a fundamental tool in time series analysis. They describe how a dynamic system reacts to a one-time shock and how this effect propagates through time.
This project illustrates how impulse responses behave under different stochastic environments using simple autoregressive processes. In particular, the analysis compares IRFs when the underlying process is stationary, unit-root (boundary case), or non-stationary.
The objective is to demonstrate how the dynamic propagation of shocks depends on the stability properties of the time series, using AR(3) processes implemented in Python.
The Stationary AR(3) Process
We consider the following autoregressive process of order three:
$ y_t = \phi_1 y_{t-1} + \phi_2 y_{t-2} + \phi_3 y_{t-3} + \varepsilon_t $
where $\varepsilon_t$ is a white noise shock with mean zero and constant variance.
In this example we use the parameters
$ \phi_1 = 1.0, \quad \phi_2 = 0.2, \quad \phi_3 = -0.3 $
The simulated time series is shown below.
Show code
import numpy as np
import matplotlib.pyplot as plt
# -----------------------------
# Parameters
# -----------------------------
phi1 = 1.0
phi2 = 0.2
phi3 = -0.3
T = 250
burn_in = 100
sigma = 1.0
seed = 42
# -----------------------------
# Simulate AR(3)
# y_t = phi1*y_{t-1} + phi2*y_{t-2} + phi3*y_{t-3} + eps_t
# -----------------------------
rng = np.random.default_rng(seed)
eps = rng.normal(0, sigma, T + burn_in)
y = np.zeros(T + burn_in)
for t in range(3, T + burn_in):
y[t] = (
phi1 * y[t-1]
+ phi2 * y[t-2]
+ phi3 * y[t-3]
+ eps[t]
)
# Drop burn-in
y_plot = y[burn_in:]
time = np.arange(1, T + 1)
# -----------------------------
# Plot
# -----------------------------
plt.figure(figsize=(8, 4.5))
plt.plot(time, y_plot, linewidth=2)
plt.axhline(0, linestyle="--", linewidth=1)
plt.title(r"Simulated AR(3) Process: $\phi_1=1.0,\ \phi_2=0.2,\ \phi_3=-0.3$", fontsize=15)
plt.xlabel("Time", fontsize=12)
plt.ylabel(r"$y_t$", fontsize=12)
plt.grid(True, alpha=0.25)
plt.tight_layout()
plt.show()
<Figure size 800x450 with 1 Axes>
Recursive Method
The first way to compute the impulse response function is by directly using the recursive definition of the AR(3) process.
We consider the model
$ y_t = \phi_1 y_{t-1} + \phi_2 y_{t-2} + \phi_3 y_{t-3} + \varepsilon_t $
To obtain the impulse response, we introduce a one-time shock at time $t=0$ and set all other shocks to zero:
$ \varepsilon_0 = 1, \quad \varepsilon_t = 0 \quad \text{for } t>0 $
Starting from zero initial conditions, the values of $y_t$ are then computed recursively using the AR equation. The resulting sequence represents the impulse response function, showing how the effect of the initial shock propagates through the system over time.
Show code
# Parameters
phi1, phi2, phi3 = 1.0, 0.2, -0.3
T = 20
# Shock at time 0
w = np.zeros(T + 1)
w[0] = 1.0
# IRF container
y = np.zeros(T + 1)
# Recursive computation
for t in range(T + 1):
y_l1 = y[t-1] if t-1 >= 0 else 0.0
y_l2 = y[t-2] if t-2 >= 0 else 0.0
y_l3 = y[t-3] if t-3 >= 0 else 0.0
y[t] = phi1 * y_l1 + phi2 * y_l2 + phi3 * y_l3 + w[t]
# Horizons
h = np.arange(T + 1)
# Plot
plt.figure(figsize=(8, 4.5))
plt.plot(h, y, marker="o", linewidth=2)
plt.axhline(0, linestyle="--", linewidth=1)
plt.title("Impulse Response Function (Recursive Method)", fontsize=14)
plt.xlabel("Horizon", fontsize=12)
plt.ylabel("Response", fontsize=12)
plt.xticks(h)
plt.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()
<Figure size 800x450 with 1 Axes>
Companion Matrix Method
The impulse response function can also be computed using the companion matrix representation of the AR(3) process.
The autoregressive model
$ y_t = \phi_1 y_{t-1} + \phi_2 y_{t-2} + \phi_3 y_{t-3} + \varepsilon_t $
can be written in state-space form as
$ x_t = F x_{t-1} + G \varepsilon_t $
where the transition matrix is
$ F = \begin{bmatrix} \phi_1 & \phi_2 & \phi_3 \ 1 & 0 & 0 \ 0 & 1 & 0 \end{bmatrix} $
and the shock loading vector is
$ G = \begin{bmatrix} 1 \ 0 \ 0 \end{bmatrix} $.
In this representation, the impulse response at horizon $j$ is obtained from
$ IRF(j) = e_1^\top F^j G $
where $e_1 = (1,0,0)^\top$.
Therefore, computing the impulse response reduces to repeatedly taking powers of the transition matrix $F$.
Show code
# Parameters
phi1 = 1.0
phi2 = 0.2
phi3 = -0.3
T = 20
# Companion matrix
F = np.array([
[phi1, phi2, phi3],
[1.0, 0.0, 0.0],
[0.0, 1.0, 0.0]
])
# Shock loading vector
G = np.array([[1.0], [0.0], [0.0]])
# Compute IRF
h = np.arange(T + 1)
irf = []
for j in h:
Fj = np.linalg.matrix_power(F, j)
irf.append((Fj @ G)[0, 0])
irf = np.array(irf)
# Plot
plt.figure(figsize=(8, 4.5))
plt.plot(h, irf, marker="o", linewidth=2)
plt.axhline(0, linestyle="--", linewidth=1)
plt.title("Impulse Response Function (Companion Matrix Method)", fontsize=14)
plt.xlabel("Horizon", fontsize=12)
plt.ylabel("Response", fontsize=12)
plt.xticks(h)
plt.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()
<Figure size 800x450 with 1 Axes>
Eigenvalue Method
A third way to compute the impulse response relies on the eigenvalues of the companion matrix.
Let $\lambda_1, \lambda_2, \lambda_3$ denote the eigenvalues of the transition matrix $F$.
The impulse response can then be written as a linear combination of powers of these eigenvalues:
$ IRF(j) = \sum_{i=1}^{3} c_i \lambda_i^{,j} $
where the coefficients $c_i$ depend on the eigenvalues and are given by
$ c_i = \frac{\lambda_i^{,p-1}} {\prod_{k \neq i} (\lambda_i - \lambda_k)} $
with $p=3$ for the AR(3) process.
In this representation, the entire dynamic behavior of the system is determined by the eigenvalues of the companion matrix. Their magnitude governs how quickly the impulse response decays, while complex eigenvalues generate oscillatory patterns in the response.
Show code
# Parameters
phi1 = 1.0
phi2 = 0.2
phi3 = -0.3
T = 20
# Companion matrix
F = np.array([
[phi1, phi2, phi3],
[1.0, 0.0, 0.0],
[0.0, 1.0, 0.0]
])
# Eigenvalues
lam = np.linalg.eigvals(F)
print("Eigenvalues of the companion matrix:")
for i, val in enumerate(lam, start=1):
print(f"λ{i} = {val:.6f}")
p = 3
# Compute coefficients c_i
c = []
for i in range(p):
num = lam[i]**(p-1)
den = 1
for k in range(p):
if k != i:
den *= (lam[i] - lam[k])
c.append(num / den)
c = np.array(c)
# Compute IRF
h = np.arange(T + 1)
irf = []
for j in h:
irf_j = np.sum(c * lam**j)
irf.append(irf_j.real)
irf = np.array(irf)
# Plot
plt.figure(figsize=(8,4.5))
plt.plot(h, irf, marker="o", linewidth=2)
plt.axhline(0, linestyle="--", linewidth=1)
plt.title("Impulse Response Function (Eigenvalue Method)", fontsize=14)
plt.xlabel("Horizon")
plt.ylabel("Response")
plt.xticks(h)
plt.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()
Eigenvalues of the companion matrix: λ1 = -0.515728+0.000000j λ2 = 0.757864+0.085703j λ3 = 0.757864-0.085703j
<Figure size 800x450 with 1 Axes>
Short summary about the comparison of methods
All three approaches — recursive simulation, companion matrix representation, and the eigenvalue method — produce identical impulse response functions.
This confirms that, despite their different formulations, the methods are mathematically equivalent and describe the same underlying dynamics of the system.
The Unit-Root AR(3) Process
In this example, we construct an AR(3) process by specifying its eigenvalues directly:
$ \lambda_1 = 1, \quad \lambda_2 = 0.3, \quad \lambda_3 = -0.8 $
These eigenvalues define the dynamic properties of the system. In particular, the presence of $\lambda = 1$ implies a unit root, meaning that the process is not stationary.
From these eigenvalues, the corresponding AR(3) coefficients are obtained as
$ \phi_1 = 0.500, \quad \phi_2 = 0.740, \quad \phi_3 = -0.240 $
Using these parameters, we simulate the time series and plot its evolution. Although the series may appear stable over short horizons, the unit root implies that shocks have permanent effects, and the process does not revert to a fixed mean over time.
Show code
# -----------------------------
# Given eigenvalues
# -----------------------------
lam = np.array([1.0, 0.3, -0.8])
# -----------------------------
# Compute AR(3) coefficients
# (from polynomial expansion)
# -----------------------------
phi1 = np.sum(lam)
phi2 = -(lam[0]*lam[1] + lam[0]*lam[2] + lam[1]*lam[2])
phi3 = lam[0]*lam[1]*lam[2]
print("AR(3) coefficients:")
print(f"phi1 = {phi1:.3f}")
print(f"phi2 = {phi2:.3f}")
print(f"phi3 = {phi3:.3f}")
# -----------------------------
# Simulate AR(3)
# -----------------------------
T = 250
burn_in = 100
sigma = 1.0
seed = 42
rng = np.random.default_rng(seed)
eps = rng.normal(0, sigma, T + burn_in)
y = np.zeros(T + burn_in)
for t in range(3, T + burn_in):
y[t] = (
phi1 * y[t-1]
+ phi2 * y[t-2]
+ phi3 * y[t-3]
+ eps[t]
)
# Drop burn-in
y_plot = y[burn_in:]
time = np.arange(1, T + 1)
# -----------------------------
# Plot
# -----------------------------
plt.figure(figsize=(8, 4.5))
plt.plot(time, y_plot, linewidth=2)
plt.axhline(0, linestyle="--", linewidth=1)
plt.title("AR(3) Process with Unit Root (λ = 1, 0.3, -0.8)", fontsize=14)
plt.xlabel("Time")
plt.ylabel("y_t")
plt.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()
AR(3) coefficients: phi1 = 0.500 phi2 = 0.740 phi3 = -0.240
<Figure size 800x450 with 1 Axes>
Impulse Response of the Unit-Root AR(3) Process
We now compute the impulse response function of the unit-root AR(3) process using the eigenvalue method.
Since one eigenvalue is equal to $1$, the response to a one-time shock does not die out completely. This reflects the fact that shocks have a permanent component in the unit-root case.
Show code
# Given eigenvalues
lam = np.array([1.0, 0.3, -0.8], dtype=complex)
T = 20
p = 3
# Compute coefficients c_i
c = []
for i in range(p):
num = lam[i]**(p-1)
den = 1
for k in range(p):
if k != i:
den *= (lam[i] - lam[k])
c.append(num / den)
c = np.array(c)
# Compute IRF
h = np.arange(T + 1)
irf = []
for j in h:
irf_j = np.sum(c * lam**j)
irf.append(irf_j.real)
irf = np.array(irf)
# Plot
plt.figure(figsize=(8, 4.5))
plt.plot(h, irf, marker="o", linewidth=2)
plt.axhline(0, linestyle="--", linewidth=1)
plt.title("Impulse Response Function (Eigenvalue Method)", fontsize=14)
plt.xlabel("Horizon")
plt.ylabel("Response")
plt.xticks(h)
plt.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()
<Figure size 800x450 with 1 Axes>
The Non-Stationary AR(3) Process
We construct an AR(3) process using the eigenvalues
$ \lambda_1 = 1.5, \quad \lambda_2 = -0.8, \quad \lambda_3 = -0.5 $
Since one eigenvalue satisfies $|\lambda| > 1$, the process is non-stationary and explosive.
The corresponding coefficients are
$ \phi_1 = 0.200, \quad \phi_2 = 1.550, \quad \phi_3 = 0.600 $
The simulated time series is shown below. Due to the explosive root, the series grows rapidly in magnitude, so we focus on a short time horizon to make the dynamics visible.
Show code
import numpy as np
import matplotlib.pyplot as plt
# -----------------------------
# Given eigenvalues
# -----------------------------
lam = np.array([1.5, -0.8, -0.5])
# -----------------------------
# Compute AR(3) coefficients
# -----------------------------
phi1 = np.sum(lam)
phi2 = -(lam[0]*lam[1] + lam[0]*lam[2] + lam[1]*lam[2])
phi3 = lam[0]*lam[1]*lam[2]
print("AR(3) coefficients:")
print(f"phi1 = {phi1:.3f}")
print(f"phi2 = {phi2:.3f}")
print(f"phi3 = {phi3:.3f}")
# -----------------------------
# Simulate AR(3)
# -----------------------------
T = 15
sigma = 1.0
seed = 42
rng = np.random.default_rng(seed)
eps = rng.normal(0, sigma, T)
y = np.zeros(T)
for t in range(3, T):
y[t] = (
phi1 * y[t-1]
+ phi2 * y[t-2]
+ phi3 * y[t-3]
+ eps[t]
)
# -----------------------------
# Plot
# -----------------------------
time = np.arange(T)
plt.figure(figsize=(8, 4.5))
plt.plot(time, y, linewidth=2)
plt.axhline(0, linestyle="--", linewidth=1)
plt.title("AR(3) Process with an Explosive Root ($\\lambda = 1.5, -0.8, -0.5$)", fontsize=14)
plt.xlabel("Time")
plt.ylabel("$y_t$")
plt.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()
AR(3) coefficients: phi1 = 0.200 phi2 = 1.550 phi3 = 0.600
<Figure size 800x450 with 1 Axes>
Impulse responses of a Non-Stationary AR(3) Process
The impulse response grows over time because the process contains an explosive root ($|\lambda| > 1$).
In this case, $\lambda_1 = 1.5$ dominates the dynamics, so the response behaves approximately like
$ IRF(j) \sim (1.5)^j $
As a result, the effect of a one-time shock not only persists, but increases over time, leading to a diverging impulse response.
Show code
# Given eigenvalues
lam = np.array([1.5, -0.8, -0.5], dtype=complex)
T = 20
p = 3
# Compute coefficients c_i
c = []
for i in range(p):
num = lam[i]**(p - 1)
den = 1
for k in range(p):
if k != i:
den *= (lam[i] - lam[k])
c.append(num / den)
c = np.array(c)
# Compute IRF
h = np.arange(T + 1)
irf = []
for j in h:
irf_j = np.sum(c * lam**j)
irf.append(irf_j.real)
irf = np.array(irf)
# Plot
plt.figure(figsize=(8, 4.5))
plt.plot(h, irf, marker="o", linewidth=2)
plt.axhline(0, linestyle="--", linewidth=1)
plt.title("Impulse Response Function of the Non-Stationary AR(3) Process", fontsize=14)
plt.xlabel("Horizon")
plt.ylabel("Response")
plt.xticks(h)
plt.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()
<Figure size 800x450 with 1 Axes>
Summary
In this project, we analyzed impulse response functions for AR(3) processes under three different regimes: stationary, unit-root, and non-stationary.
We computed IRFs using three equivalent methods: recursive simulation, companion matrix representation, and the eigenvalue approach. The results showed that all methods produce identical impulse responses.
Finally, we demonstrated how the behavior of IRFs depends on the underlying dynamics: decaying in the stationary case, persistent under a unit root, and explosive when the system is unstable.