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

Steady State Variables

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()