Let's import some scientific libraries
import numpy as np
import matplotlib.pyplot as plt
The next instruction tells Jupyter to show plots inside the browser
%matplotlib inline
Now here's our code. It plots the time series
$$ X_{t+1} = \alpha X_t + W_{t+1} $$where
and $\alpha$ takes three different values
alphas = [0.0, 0.8, 0.98]
ts_length = 200
fig, ax = plt.subplots()
for α in alphas:
X = np.zeros(ts_length) # Allocate an array of zeros
for t in range(ts_length - 1):
X[t+1] = α * X[t] + np.random.randn()
ax.plot(X, label=r'$\alpha = {}$'.format(α))
ax.legend(loc='lower left', frameon=False, ncol=3)
<matplotlib.legend.Legend at 0x7f394924cfd0>
Generate a single time series following an AR(1) process as above when $\alpha = 0.99$. However, let's now add reflecting barriers at $-2$ and $2$. Thus,
$$ y_{t+1} = \alpha X_t + W_{t+1} $$and
$$ X_{t+1} = \begin{cases} -2 & \text{ if } y_{t+1} < -2 \\ 2 & \text{ if } y_{t+1} > 2 \\ y_{t+1} & \text{ otherwise} \end{cases} $$Plot the time series.
for i in range(40):
print("solution below!")
solution below! solution below! solution below! solution below! solution below! solution below! solution below! solution below! solution below! solution below! solution below! solution below! solution below! solution below! solution below! solution below! solution below! solution below! solution below! solution below! solution below! solution below! solution below! solution below! solution below! solution below! solution below! solution below! solution below! solution below! solution below! solution below! solution below! solution below! solution below! solution below! solution below! solution below! solution below! solution below!
α = 0.99
X = np.zeros(ts_length)
for t in range(ts_length-1):
y = α * X[t] + np.random.randn()
if y > 2:
X[t+1] = 2
elif y < -2:
X[t+1] = -2
else:
X[t+1] = y
fig, ax = plt.subplots(figsize=(9, 6))
ax.set_ylim(-3, 3)
ax.plot(X, label="AR1 with reflecting barriers")
ax.legend(frameon=False)
<matplotlib.legend.Legend at 0x7f3948c47b00>
Assuming that the spectral radius of the $n \times n$ matrix $\Pi$ is strictly less than one, the equation
$$ x = \Pi x + b $$has a unique $n \times 1$ solution vector $x$ for each choice of the vector $b$. The solution satisfies
$$ x^* = (I - \Pi)^{-1} b = \sum_{t=0}^\infty \Pi^t b $$where $I$ is the identity matrix. Try computing this solution both ways, for the matrix
$$ \Pi = \begin{pmatrix} 0.4 & 0.2 \\ 0.6 & 0.2 \end{pmatrix} $$and vector $b = \begin{pmatrix} 1 \\ 1 \end{pmatrix}$. In the second case you will have to truncate the infinite sum at some reasonable level.
for i in range(40):
print("solution below!")
solution below! solution below! solution below! solution below! solution below! solution below! solution below! solution below! solution below! solution below! solution below! solution below! solution below! solution below! solution below! solution below! solution below! solution below! solution below! solution below! solution below! solution below! solution below! solution below! solution below! solution below! solution below! solution below! solution below! solution below! solution below! solution below! solution below! solution below! solution below! solution below! solution below! solution below! solution below! solution below!
First let's solve it by matrix inversion.
from numpy.linalg import solve
Π = np.array([[0.4, 0.2],
[0.6, 0.2]])
b = np.ones(2).reshape(2, 1)
x_star = solve(np.eye(2) - Π, b)
print(x_star)
[[ 2.77777778] [ 3.33333333]]
Now let's try using the sum and see if we get the same answer:
E = np.eye(2)
x_star = np.zeros((2, 1))
for i in range(100):
x_star += E @ b
E = E @ Π
print(x_star)
[[ 2.77777778] [ 3.33333333]]