In [1]:
import numpy as np
import pandas as pd
import seaborn as sns
import matplotlib.pyplot as plt

import scipy.linalg

sns.set_style("whitegrid")


### Coefficients¶

In [2]:
y = .5  # persistence of technology LOM
g = .1
beta = .99
alpha = 0
gamma = 0


In the following variable names: 'ss' ~ 'steady state'.

In [3]:
ss_return = 1 / beta
ss_consum = y - g

rho_1 = (-beta**2 * ss_consum * alpha) / (1-beta)**2 + ss_consum / (1-beta) + 1 / beta
rho_2 = (ss_consum * alpha) / (ss_return-1)**2 - alpha
rho_3 = -ss_consum * beta**2 / (1-beta)**2
rho_4 = ss_consum * beta**2 / (1-beta)**2 - 1

tmp0 = (gamma - 1 / beta)
tmp1 = tmp0 * rho_3 - rho_4
tmp2 = tmp0 * rho_1 - rho_2


### Dynamic System¶

The matrices can be constructed in a similar fashion as on the last problem set.

In [4]:
A = np.array([
[1, rho_1, 0, 0],
[0, 1    , 0, 0],
[0, 0    , 1, 0],
[0, 0    , 0, 1],
])

B = np.array([
[-tmp0, -rho_2    , tmp1 , tmp2      ],
[0    , alpha*beta, beta , alpha*beta],
[0    , 0         , 0    , 0         ],
[0    , 0         , 0    , 0         ],
])

C = np.array([tmp0, 0, 0, 0])


In the function scipy.linalg.ordqz I set sort="ouc" which stands for outside unit circle. This achieves the desired order for our problem. (Here I am reusing the code written for the 4th problem set; see here.)

In [5]:
S, T, _, _, Q, Z =  scipy.linalg.ordqz(A, B, sort="ouc")

In [6]:
eigen_values = np.diag(T) / np.diag(S)
eigen_values = np.round(eigen_values, decimals=15)

I = np.where(np.abs(eigen_values) < 1)[0]

ns = len(I)

S_11, T_11, Z_11 = (m[:ns, :ns] for m in (S, T, Z))
Z_21 = Z[ns:, :ns]

Z_11_inverse = np.linalg.inv(Z_11)

In [7]:
S2S = Z_11 @ np.linalg.inv(S_11) @ T_11 @ np.linalg.inv(Z_11) # states to states
S2C = Z_21 @ np.linalg.inv(Z_11)                              # states to controls


### Actual Computation of the Impulse Response¶

In [61]:
n_periods = 5
state = np.array([0, 0, 1])  # 1 corresponds to 1% shock
result = np.empty((n_periods, 8))

In [62]:
for t in range(n_periods):
controls = S2C @ state
updated_state = S2S @ state
nominal_rate = beta * state[2]

result[t, :ns] = state
result[t, ns] = controls
result[t, ns + 1] = controls[0] + state[1]  # inflation
result[t, ns + 2] = state[0] - rho_1 * controls[0] - rho_3 * state[2]  # bonds
result[t, ns + 3] = nominal_rate
result[t, ns + 4] = rho_3 * nominal_rate  # money demand

state = updated_state

result = np.round(result, decimals=12)  # remove numerical annoyances


### Plotting¶

In [65]:
df = pd.DataFrame(
result[:, 2:],
columns=[
"monetary_shock",
"surprise_inflation",
"inflation",
"bonds",
"nominal_rate",
"money_demand"
],
index=np.arange(n_periods)
)
df.index.name = "periods"

In [66]:
for col in df.columns:
fig, ax = plt.subplots(1)
fig.set_size_inches(15, 8)

sns.lineplot(
x="periods",
y=col,
data=df[[col]],
linewidth=2,
ax=ax,
)

ax.set(xlabel="Periods")
ax.set_frame_on(False)
ax.xaxis.label.set_size(13)
ax.yaxis.label.set_size(13)
plt.show()